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1. Introduction 


The mathematical formulation of the engineering optimization problem is 

min /*(*}) 

subject to g,({x})i;0, i=l,q 


( 1 ) 


where 

{x} is an nxl matrix of design variables, 
f({x}) is the objective function, and 
gj({x}) are constraint equations. 

Evaluation of the objective function and constraint equations in Equation (1) can be very 
expensive in a computational sense. Thus, it is desirable to use as few evaluations as 
possible in obtaining its solution. In solving Equation (1), one approach is to develop 
approximations to the objective function and/or restraint equations and then to solve 
Equation (1) using these approximations in place of the original functions. These 
approximations are referred to as response surfaces. 

The desirability of using response surfaces depends upon the number of functional 
evaluations required to build the response surfaces compared to the number required in the 
direct solution of Equation (1) without approximations. The present study is concerned with 
evaluating the performance of response surfaces so that a decision can be made as to their 
effectiveness in optimization applications. In particular, this study focuses on how the 
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quality of approximations is effected by design selection. Polynomial approximations and 
neural net approximations are considered. 

To provide the groundwork for future discussion, this introductory section discusses: 

1. measures of quality of fit at the designs and measures of quality of fit over a region of 
interest and 

2. the methodology used to build the approximations. 



Let us consider a problem with n design variables, the components of the vector {x} = {x 1? 
x 2 ,...x n }‘. A total of N designs will be considered: {x} jt j = l,N. At the designs {x } j5 let 


y. = the value of the function to be approximated and 
yj = the value of the approximating function. 

The approximating function, % should closely match the function, y, not only at the designs, 


{x}j, but over the entire region of interest. 



The approximating function $ closely approximates the function y when s is small where 



and where S 2 is the sum of the squares of the residuals thus 


2 


( 3 ) 


6 2 =T(y r Sf 

1 

Let y be the average value of the designs, y,. Thus 


N 



(4) 


In this study, one measure of the closeness of fit to be considered is the non-dimensional 
value v where 


\ 

V= — 




N 


* 100 


(5) 


The coefficient v is the non-dimensional root mean square (RMS) error at the designs. 
Thus, v = 0 is a necessary and sufficient condition that the approximating function fit the 
actual function at the N design points. 

1.1.2 Overall fit 

Just because the approximating function exactly fits the function at N designs does not 
guarantee that it gives a good fit over the region of interest. It is therefore desirable over 
the region of interest to have a measure of the quality of overall fit. Several examples of 
this study considers a two dimensional region of interest. For these problems, the 
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rectangular region of interest is overlaid with a 31x31 evenly spaced grid of points. The 
value of the function and the approximating function is then compared at these NG = 961 
evenly spaced grid of points. Other examples consider a rectangular n dimensional region 
of interest. These regions of interest are also overlaid with a evenly spaced grid of points. 
The value of the function and the approximating function are then compared at these NG 
grid points. For these examples, a measure of the quality of overall fit is taken as 


v c=- 


NG 

^ ZiyrSflNG 

Tg 


* 100 


( 6 ) 


where y G is the average value of y at the grid points. A small value of v G indicates that the 
approximating function did a good job of approximation over the region of interest. 

1.2. Polynomial Approximations 

With the polynomial response surface approach, the approximating function is taken as an 
m=k + 1 term polynomial expression [1-3] thus 

y=b a +b 1 X 1 +...b i X k (7) 

where Xj is some expression involving the design variables. For example, a second order 
polynomial approximation in two variables could be of the form 
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y= b o +b l X l +b 2 t 2 +b 3 X l +b ^ C l X 2 +b sh (8) 

The value of the function to be approximated at the N designs can be used to determine the 
m = k+ 1 undetermined coefficients in the polynomial expression. For the N designs, 
Equation (7) yields 




[1 X, ... X. 

K 

y 2 


1 X . ... X. 

l i 


< 

y 



y* 


1 X, ... x t 

K 


(9) 


or 

1 Y\=[Z\\b\ ( 10 > 

where {Y} is an Nxl matrix, [Z] is an Nxm matrix, and {b} is an mxl matrix. 

1.2.1 Exactly-determined approximation 

When N = m, the approximation is exactly-determined and the matrix {b} can be determined 
from Equation (10). 

1.2.2 Over-determined approximation 

With N>m, Equation (10) can be solved in a least squares sense thus [1-3] 


5 



izfm-mznb) 


(id 


or 

w-aznzD-wm (12) 

Equation (12) in effect, chooses the terms of {b} so as to minimize the square of the 
residual as defined in Equation (2). 

1.2.3 Under-determined approximation 

When N<m, the approximation is under-determined. A solution can be obtained by 
choosing the terms of {b} so as to minimize the square of the residual as defined in 
Equation (2). However, a direct solution can be obtained by using the concept of pseudo- 
inverse [4,5]. Assume that the rank of matrix [Z] is N and define the pseudo-inverse of 
matrix Z, Z’ thus 

[zjMZMZHzrr 1 (13) 

where t denotes transpose. Solution of Equation (10) is then 

{M=[Zrm+K?]{W (14) 
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where {w} is an (m-N) column matrix of arbitrary coefficients and [Q] is a mx(m-N) matrix 
formed from any m-N independent columns of the matrix [R] thus 

[K\=[I]-[Z\'[Z\ < 15 > 

One solution to Equation (14) is to take all the arbitrary terms of (w) as zero giving 

{frMznn (i6) 

The basic solution to Equation (10) is Equation (16). Using that equation, at the designs, 
{x}j, the value of matches the value of y jt If w ( is the ith term in matrix {w} and {q}; is 
the ith column of matrix [Q], then at the designs, {x} jf ^=0 when 

{*}=*,{$ 1, (17) 

Thus, the last term of the right hand side of Equation (14) gives ^ values which match 
at the designs, {x}j, for any values of Wj. 

1.3 Artificial Neural Nets 

While the initial motivation for developing artificial neural nets was to develop computer 
models that could imitate certain brain functions, neural nets can be thought of as another 
way of developing a response surface. Different types of neural nets are available [6,7], but 
the type of neural nets considered in this paper are back propagation nets with one hidden 
layer as shown in Figure 1. This type of neural net has been used previously to develop 
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response surfaces [8-12] and is capable, with enough nodes on the hidden layer, of 


appro ximating any continuous function [13]. 


For the neural net of Figure 1, associated with each node on the hidden layer, node j, and 
each output node, node k, are coefficients or weights, 0j and 0^ respectively. These weights 
are referred to as the biases. Associated with each path, from an input node i to node j on 
the hidden layer, is an associated weight, w^ and from node j on the hidden layer to output 
node k is an associated weight w jlc . Let q be inputs entered at node i. Node j on the 
hidden layer receives weighted inputs, w^. It sums these inputs and uses an activation 
function to yield an output r^ The activation function considered in this paper is the 
sigmoid function [6,7] 



J 

-Ew # fl r e / 


( 18 ) 


Output node k then receives inputs which are summed and used with an activation 
function to yield an output s k . Some variation of the delta-error back propagation algorithm 
[6,7] is then used to adjust the weights on each learning try so as to reduce the values 
between the predicted and desired outputs. In this investigation, studies were performed 
using the program NEWNET [14] which was developed especially for this investigation. 
NEWNET minimizes the sum of the squares of the residuals in Equation (2) with respect 
to the weights and biases of the net. Training of the net is thus formulated as an 
unconstrained minimization problem. Solution of this minimization problem is performed 
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using the method of Davidon, Fletcher, and Powell [15-16]. That algorithm performs a 
series of one dimensional searches along search directions. Search directions are 
determined by building an approximation to the inverse Hessian matrix using gradient 
information. Gradients required by that algorithm are obtained using back-propagation. 
One-dimensional searches are performed along the search directions using an interval 
shortening routine. 
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2. Levels of Designs 


2.1 Tavlor Series Approximation 

The overriding factor which affects the accuracy of an approximation is the levels of the 
design parameters considered. It is instructive to consider a problem in two design 
variables. Suppose we wish to make a quadratic approximation of a function thus: 

y=b 0 +b 1 x l +bjX 2 +b 3 x l +b 4 x l x 2 +b s x 2 ... ^ ^ 


Consider that the exact function is evaluated at 6 design points and the information thus 
generated will be used to determine the 6 undetermined coefficients in Equation (19). 
Design variables at these design points are taken from the following sets. 

Xj from the set x^ (20) 
x 2 from the set {x 2j x^-.-x^} 

Here p discrete values are considered for Xj and q discrete values are considered for x 2 . 
The variable is said to have p levels and x 2 is said to have q levels. The problem is to 
determine the minimum levels of the design variables, p and q, required to build the 
quadratic approximation. In this regard, it is instructive to consider a Taylor series 
approximation [17] of the function about the point {xj=0, x 2 =0}: 

y=y(0,0)+ 1 vy(0,0)}'{Ax} + {Ax} r [tf(0,0)]{/Lx} + . . . (21) 
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where 


{Ax}=[(jc 1 -o)(x 2 -o)r=[jt 1 x 2 r 


( 22 ) 


{vy(0,0)}=[( 


M °>°) MMy 

dx x 


( 23 ) 


[^(0,0)]= 


a^Q.O) ay 2 (0,0) 

dxi 

ay 2 (Q,Q) ^y(Q,Q) 


dXjCbtj SXj 


( 24 ) 


Entering Equations (22), (23), and (24) into Equation (21) gives 


y=y(0fi) 



( 25 ) 


The derivatives in Equation (25) can be determined by finite difference equations [18]. The 
second derivative of y with respect to x, can be obtained using information at points 
indicated in Figure 2 by solid circles, the second derivative of y with respect to x 2 can be 
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obtained using information at points indicated by unfilled circles, and the mixed derivative 
can be obtained using information at points indicated by unfilled squares. 

It can be seen in Figure 2 that at least three levels of both X! and x 2 must be used to obtain 
a quadratic approximation. If three levels are not provided, not information is available to 
calculate the higher derivatives in Equation (25). A complete 3 factorial design does not 
have to be used-only 6 selected points from the complete 3 factorial design. Information 
at those 6 points allow the undetermined coefficients to be exactly determined. 

Consider now the design of Figure 3 which are also taken from the 3 factorial design. Even 
though 6 design points are used, this set of design points does not allow an approximation 
containing the x 2 2 term of Equation (25). However, with the design of Figure 3, an 
approximation of the form of Equation (26) could be obtained thus: 

y=VVi + V2 + Vi 2+ *w 2 (26) 

With the design of Figure 3, if a solution is attempted using Equations (19) and (12), a 
singular coefficient matrix will be encountered. A solution could be attempted using the 
pseudo-inverse concept of Equations (13) and (14). However, recent studies [19] have 
shown that non-unique solutions are obtained with this technique. Non-uniqueness makes 
these solutions undesirable. Using Equations (26) and (12), a slightly over-determined 
approximation is obtained. 
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Recent studies have found that the numerical performance of neural network 
approximations and polynomial approximations with the same number of associated 
undetermined parameters is comparable [19]. Thus, it is not expected that neural nets as 
approximators will perform better than polynomials when there are inadequacies in the 
training design, as in Figure 3. The next example investigates performance of both 
polynomial and neural net approximations. 

2.2 Example 

Consider the function 


y = 1 +x i +x 2 +x 3 +X 1 + X 1 X 2 +X 1 X J +X? -t-XjXj+Xj 2 


(27) 


In the first phase of the investigation, approximations are to be made of this function using 
the design of Figure 4. The star pattern of design points in Figure 4 does not allow mixed 
derivatives of the function to be calculated using finite difference type formulae but does 
permit the other second derivatives to be calculated. Thus, information is available to make 
a polynomial approximation of the form 

y=b o +b l x l +b 2 x 2 +b 3 x 3 +b 4 x* +b s x%+b 6 x* (28) 

The function y was evaluated at the design points shown in Figure 4 yielding 7 training pairs 
for calculating the 7 undetermined parameters in Equation (28). The value of the 
approximating function y was then evaluated at a 5x5x5 grid of designs. These values of $ 
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were then used to evaluate v G from Equation (6). The value of v G obtained is shown in the 


first line of Table 2.1. 


Table 2.1. Performance of Approximations for Various Designs 


Number 

Designs 

Points 

Description 

Polynomial 

Approximation 

Neural Net 
Approximation 

No. 

Para. 

v G (%) 

ih 

No. 

Para. 

No. 

Apx. 

v G (%) 

7 

Star-see 
Figure 4 

7 

34.6 

2 

11 

10 

25.5-97.3 

12 

Star-see 
Figure 5 

7 

34.6 

2 

11 

10 

32.9-93.5 

10 

Computer 

Generated 

10 

0.0 

2 

11 

10 

36.6-36.9 

3 

16 

10 

21.9-36.7 

27 

3 factorial 

10 

0.0 

3 

16 

2 

16.6-16.7 

a 

21 

2 

16.6-16.9 

125 

5 factorial 

10 

0.0 

8 

41 

1 

3.7 


A neural net approximation was then considered. Previous studies [19] have indicated that 
it is desirable to have more training pairs than the number of undetermined parameters 
(weights and biases) associated with the net. If fewer training pairs than undetermined 
parameters are used, non-unique approximations should be expected. For a neural net with 
one hidden layer as shown in Figure 1, there are 6 parameters associated with a net with 
one node on the hidden layer and 11 parameters associated with a net with two nodes on 
the hidden layer. It was considered that one node on the hidden layer would yield an 
inadequate approximation. Thus 2 nodes on the hidden layer were considered. Thus, the 
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neural net approximation is under-determined. That is to say that there are fewer training 
pairs than there are undetermined parameters associated with the approximation. Non- 
unique approximations are to be expected. Indeed, this was the case. The 8 training pairs 
were used to make 10 different approximations by having training commence from a 
different randomly selected set of weights and biases. Once the nets were trained, the value 
of the approximating function, % was generated at the 5x5x5 set one grid points and the 
value of v G was developed. The range of the values obtained is shown in Table 2.1. One 
can see that a large range of values is obtained. The best neural net approximation is only 
slightly better than the polynomial approximation while the worst neural net approximation 
is considerably worse. Just as with the polynomial approximation, the designs used to train 
the approximation can not yield information necessary to capture essential features of the 
function to be approximated. 

The 12 designs of Figure 5 were next used in the training of a polynomial approximation 
and a 2 node neural net approximation. Even though more designs are used here than in 
Figure 4, the additional designs selected do not yield any more information about the nature 
of the function being approximated. Information is still not available for determining the 
mixed derivatives of the function to be approximated. Thus, the polynomial approximation 
of Equation (26) was considered. As there are now more training pairs than there are 
undetermined parameters, the approximation obtained is over-determined. As no new 
information is available with the 12 designs, the same polynomial approximation and thus 
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the same v G as before are obtained. The value of v G is shown in the second line of Table 

2 . 1 . 


A neural net with 2 nodes on the hidden layer was then trained with the 12 training pairs. 
The net was trained 10 times starting from different randomly selected sets of weights and 
biases. Even thought the number of training pairs, 12, is greater than the number of 
undetermined parameters associated with the net, 11, non-unique approximations were 
obtained as can be seen in Table 2.1. Thus, it can be concluded that for neural net 
approximations, having more training pairs than the number of associated undetermined 
parameters is only a necessary condition for obtaining a unique approximation but that it 
is not a sufficient condition . As the 12 designs offered no new information about the 
function being approximated over that offered by the 8 designs, then just as with the 8 
design case, non-unique approximations were obtained. 

The program DESIGNS [20], which was developed for this project, was used to generate 10 
designs which contain the information necessary for calculating the 10 undetermined 
coefficients of the complete quadratic approximation of the form: 

y= b o +b l X l + V 2 +fe 3*3 + Vl 2 + V2 + Vs 2 + VA + V 1*3 + V2*3 (29) 

The location of these design points is shown in Figure 6. The polynomial approximation 
obtained by training the polynomial of Equation (29) with the computer generated designs 
exactly duplicated the test function of Equation (27). Thus, v G for the 5x5x5 grid of points 
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was zero as seen in the third line of Table 2.1. 


A neural net with 2 nodes on the hidden layer with 6 associated undetermined parameters 
and a neural net with 3 nodes on the hidden layer and 11 associated undetermined 
parameters were then trained 10 times with the computer generated training pairs. Each 
training started from a different randomly selected set of weights and biases. For the case 
of 2 nodes on the hidden layer, the approximation generated was over-determined and a 
unique approximation was obtained (the small range of v G obtained most likely results from 
the exit criteria employed in the training algorithm). For the case of 3 nodes on the hidden 
layer, there are 11 associated undetermined parameters but only 10 training pairs. Thus the 
approximation is under-determined and a non unique approximation is obtained as can be 
seen in Table 2.1. 

The performance of the neural net approximations was much poorer than that of the 
polynomial approximation on this problem. This poorer performance may be in part 
because the problem is biased towards the polynomial approximation as the function being 
approximated is 2 second order polynomial. 

A complete 3 3 factorial design and a 5 3 factorial design were considered to see if good 
results could be obtained with the neural nets if more training pairs were employed. Indeed 
this was the case. However, many more training pairs were required to get a good 
approximation than were required with the polynomial approximation. The extra training 
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pairs were wasted on the polynomial approximation. Ten correctly selected training pairs 
is all that is required to get an exact second order approximation. The additional training 
pairs offered no new information to the polynomial approximation. The coefficient v G was 
zero for training pairs using the 3 and 5 factorial designs and a second order polynomial 
approximation. 

2.3 Conclusion 

For a given order of approximation, a good design must use an adequate number of levels 
of the design variables or a poor approximation will be obtained. Likewise, design points 
must be located so that information is available for determining all of the undetermined 
coefficients of the approximating function. In many instances, especially when the region 
of interest is small, a second order polynomial approximation or neural net equivalent will 
be sufficient to build a response surface. A second order approximation requires a design 
containing 3 levels of the design variables. Program DESIGNS has been developed to 
generate a minimum point design which allows all of the coefficients of a second order 
polynomial approximating function to be obtained. This minimum point design can be 
augmented by randomly selected design points or by user selected points. 
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3. Standard Designs 


3.1 Underlying Principle 

When making a polynomial approximation of a function, the number of design levels 
required for each design variable depends upon the order of polynomial approximation 
being used. Consider for example the problem of approximating a function y, a function of 
one design variable. As previously discussed, two levels of the design variable would be 
required to make a linear approximation of the function, three levels of the design variable 
would be required to make a second order approximation, four levels of the design variable 
would be required to make a 3rd order approximation, etc. If y is a function of r design 
variables, a pth order polynomial approximation , % requires designs at p+1 levels in each 
design variable. 

In response surface methodology, the term factor is used for design variable. A factorial 
design or factorial experiment is a design in which one uses each of the possible 
combinations of the levels of each factor. If m is the number of level of each factor and r 
is the number of factors, then the design would be referred to as a m r factorial experiment . 
Table 3.1 gives the number of designs in various factorial experiments. 
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Table 3.1. Number of designs in a full factorial design 


m = level 
r = factor 

2 

3 

4 

2 

4 

9 

16 

3 

8 

27 

64 

4 

16 

81 

256 

10 

1024 

59049 

1.05E06 


One can see that even for a small number of factors, complete factorial experiments become 
impractical if designs are computationally or experimentally expensive to obtain. One then 
is forced to use some sub-set of the factorial design or alternate designs containing requiring 
fewer design points. Concepts from statistics are normally used in selecting a sub-set of the 
factorial design or in developing alternate designs. Thus statistical concepts are reviewed. 

3.2 Statistical Concepts 

When making an approximation, ?, of a function, y, most approaches used to select design 
points for a design consider that 

1. polynomial approximations are employed and 

2. the value of the function, y,, determined at the designs, {x} ; , contains some error, e { . 
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A measure of the error at point i is the variance of the error, var(€j) = a 2 where 


A 



i« 1 


(y,-^ z 


( 30 ) 


where 

M is the true mean of all possible observations of y ; and 
n is the number of observations made. 

In experimental investigations, is experimental error. When making approximations to 
analytical functions, e i is zero and the variance of the error at point i is zero. Often 
approximations are made to a function whose values must be obtained from some numerical 
algorithm such as the finite element method or finite difference method. Values of y ; 
obtained from such algorithms depend on control parameters which dictate the level of 
accuracy of the solution. For example, if y was a stress determined from a finite element 
analysis, then y could depend on a control parameter which specifies the coarseness of the 
finite element idealization. In this case, different values of y { would be obtained for the ith 
design for different values of the control parameters and e t could be thought of as a 
numerical error. 

It would be an interesting study to select designs such that approximations developed are 
insensitive to numerical errors such as finite element idealization error. However, the 
problem at hand is to find a good approximation to an analytical function or a good 
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approximation for output from a deterministic model. For the problem at hand, for a given 
design, Xj, one obtains the same functional value, y ( no matter how many times the function 
is evaluated. Thus, the problems considered in this report contain no numerical error. 
However, as all known algorithms with one exception [21] consider that there is some 
experimental or numerical error, this section now further examines this case. 

Errors in the value of y, used to build an approximation affect the estimation of the 
undetermined coefficients, bj, in the polynomial approximation and thus affect ft, the values 
of y ; predicted by the approximation. A measure of the error in bj resulting from errors in 
y ; is the variance of bj. For example, consider that y, is obtained from a finite element 
analysis and that a pth order polynomial approximation is employed. The undetermined 
coefficients in that approximations, bj, can be determined from Equation (12). If a number 
of approximations were now made with finite element results, obtained using different 
idealizations, the coefficient bj for these approximations would be different. The variance 
of bj is a measure of how much the b’s change for these different approximations. In like 
form, the different approximations yield different ft and the variance of ft is a measure of 
how much the ft values change from approximation to approximation. 

From a numerical standpoint, it is desirable to have approximations that are not highly 
sensitive to the error e r Approximations are insensitive to the error, e it if the variance of 
bj and the variance of ft is small. Most design selection algorithms currently in use attempt 
in some way to keep these variances small. 
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The variance of bj is the j,j term of the variance-covariance matrix cov b where (see 
Equation 3.11 of [3] or Equation 2.8 of [2]) 

[cov b] = o 2 ([Z] f [Z])'* < 31 > 

and the variance of ft is given by (see Equation 2.11 of [2]) 

vary^iZimmZir'W (32) 

where {Z,} 1 is the lxp vector whose elements correspond to the elements of a row of matrix 

[Z]. 

Notice that these variance involve the matrix [H] where 

[HMizymr 1 (33) 

Design selection affects [Z], which from Equation (33) affects [H], which in turn affects the 
variances of bj and ft. Many design point selection algorithms attempt to select designs 
which give an [H] matrix which will keep the variances of bj and ft small. 

3.3 Orthogonal Designs 

The associated undetermined coefficients of a polynomial approximation function can be 
found from Equation (12). The solution for these coefficients involve the matrix [Z] (see 
Equations (9) and (10)). Let {ZJ be the ith column of matrix [Z]. A design is said to be 
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orthogonal if the columns of the [Z] matrix are orthogonal, i.e. {ZJ'fZj} = 0, i^j. There are 
interesting properties of orthogonal designs which have prompted there use. Thus 
orthogonal designs will now to presented in some detail. 

3.3.1 Scaling 

The discussion of orthogonality is simplified by working with scaled variables. Consider that 
the approximation in question involves k unsealed design variables \ and contains N design 
points. Instead of working with x^ the variables will be scaled. Let x^ be the uth level of 
unsealed variable i and x iu be the scaled level. The desired scaling is 

i-u <30 

U m 1 


N 


E x » =0 * 


i=i jt 


This scaling can be accomplished by having 



(35) 


(36) 


where 


and 


x { =the average of the levels of x. 


(37) 
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( 38 ) 




-i N 

With this scaling, N experimental design points of the orthogonal design give 


[znz\=m 

(39) 

([Zftzjr^lm 

N 

(40) 


where [I] is the identity matrix. 


3.3. 1.1 Example of Scaled Designs: 

Consider a 2 factorial design with levels of 4 and -4. For that design 

* x =0, Xj=0 (41) 


and 


c 2 _c 2 _ (4-0) 2 +(-4-0) 2 
-*2 1 . 


or S t =S 2 = 4 


From Equation (3), the levels of the scaled variables are 




^-o 


or the levels of the scaled variables are 1 and -1. 


(42) 


(43) 
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3.3.2 Bias 

Assume that the polynomial approximating function is inadequate. The coefficients of that 
polynomial can be determined from Equation (12). Let {6j} be the coefficients thus 
obtained and let [ZJ be the corresponding [Z] matrix. Then from Equation (12) 

{^ 1 }=([Z 1 ] , [Z 1 ]) , [Z i r{n (44) 


Assume that the function being approximated can be expressed as 

{Y]=[Z]{b) 


(45) 


where 



[Z\=l [ZJ 


[ZJ 1 


(46) 


Entering Equations (40), (45), and (46) into Equation (44) gives 


(«, 


l-il/HZjW,] [Zj]) 
N 



Entering Equation (39) into Equation (47) gives 

l£j-T,WHi> 1 )*[Z 1 ]'[ZJ[fr 2 l) 

N 


(47) 


( 48 ) 


or 
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(49) 


{*, ) ={*, » +^[Z 1 ]'[Z 2 ]{fc 2 } = {b l }+[A] {b 2 } 


where [A] is called the alias matrix . One can see in Equation (49) that the coefficients {bj} 
will only be correct estimates of {bj if the columns of [ZJ are orthogonal to the columns 
of [Z 2 ], Special situations where this orthogonality occurs are next discussed. 

3.3.2. 1 A bias example-linear approximating polynomial but the exact function contains 
linear terms and cross-product terms: 

Consider a linear approximating polynomial 

< 5# > 

i«l 


where the exact function is 


H + E*A + EEVi*y 

f-1 i-l j-i 


(51) 


where b^ are the undetermined coefficients associated with the cross-product terms. For this 
problem, a full 2 k factorial design gives that the columns of [ZJ are orthogonal to the 
columns of [Z 2 ] and thus 
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( 52 ) 


33.2.2 A bias example-linear approximating function but the exact function is a complete 
quadratic polynomial: 


Consider a linear approximating polynomial 


s- 






i-1 


(53) 


where the exact function is a complete second order polynomial thus 

>*VE *E *«*?*E E Va 

i-l i-1 i-1 j-i 


(54) 


Assume again that a full 2 k factorial design is used. For this problem the alias matrix is 
such that one obtains 


K= b o+Y, b u 

i-1 

Sj=bj, j= 


Thus only 6 0 is biased with the other coefficients unbiased or uncorrelated. 
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3.3.3 Orthogonal Designs for Linear Approximations 

For a problem with r design variables, a full 2 r factorial design is an orthogonal design if the 
approximating function is a first order polynomial. There are several advantages in using 
such an orthogonal design when the approximating function is assumed to be linear. These 
advantages are: 

1. The solution for the coefficients of the polynomial approximation require a matrix 
inverse (see Equation (12)). However, when the design is an orthogonal design, that inverse 
is very easily obtained using Equation (40). Thus there is a small computational advantage 
in using an orthogonal design. 

2. Examples 3.3.2. 1 and 3.3.2.2 indicate that under certain conditions, the coefficients 
obtained using an orthogonal design are unbiased. Obtaining unbiased coefficients is 
probably more important in developing response surface from experimental results than 
when developing response surfaces when results are from a deterministic model. With 
experimental studies, it may be important to ascertain the unbiased values of the linear 
coefficients. For the deterministic model however, one is looking for an 
approximating function which gives a good approximation throughout a region of interest. 
Whether the coefficients of the polynomial approximation are biased or unbiased is of little 
concern. 

3. It can be proven that for linear polynomial approximations, an orthogonal design gives 
the minimum variance of the coefficients (see page 109 of [3]). It is important when 
modeling experimental results to obtain a model that is not overly sensitive to experimental 
error and thus there is an advantage in having a minimum variance of the coefficients. 
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However, for response surfaces of a deterministic model, variance of the coefficients is not 
relevant. 

3.3.4 Orthogonal Designs for 2nd Order Polynom ial Approximations 
It is not possible to find an orthogonal design when using a second order polynomial 
approximating function of the form of Equation (8) (see page 107 of [2]). However, an 
orthogonal design can be found if one uses as the approximating function a second order 
orthogonal polynomial (page 130 of [3]) 

i«l i»l j m i 


where 



N 



M-l 


N 


(57) 


and where 


N=the number of design points and 
x , =x, for each of the design points . 

J u J 


(58) 


The use of an orthogonal design still gives the small computational advantage that the 
inverse shown in Equation (12) is an inverse of a diagonal matrix. However, when using 
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second order approximations, it is not clear under what conditions one obtains unbiased 
coefficients. Also it can not be proven that orthogonal designs any longer give a minimum 
variance of the coefficients. Thus most of the reasons for using orthogonal designs found 
for linear approximations are not present when using second order approximations. 


H5 Ge neral Discussion of Orthogonal Designs 

Orthogonal designs offer a small computational advantage that the matrix inverse required 
m solving for the coefficients of the polynomial approximating function is an inverse of a 
diagonal matrix. When approximating a deterministic model, properties of orthogonal 
designs which minimize the variance of the coefficients and which give unbiased coefficients 
are unimportant. For this case, the use of orthogonal designs can only be justified by how 
well they perform on test problems. Such test problems are presented later in this report. 

1 4 Central Compos ite Designs -Desiens for Fitting Second Order Models 
It was shown in Section 2 that at least 3 levels of the design variables are required if one 
is to make a second order approximation. A workable alternative to using a 3 k factorial 
design is a class of designs called the central composite design . These types of designs are 
widely used by workers applying second order response surface techniques [3]. 


3,4,1 Fo rmat of the central composite design 

The central composite design is a design composed of the 2 k factorial design augmented by 
additional points. The augmented design points are as follows: 
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*1 *2 *3 "* X k 

0 0 0 ... 0 

-a 0 0 ... 0 

a 0 0 ... 0 

0 -a 0 ... 0 < 5 ’> 

0 a 0 ... 0 

»•« ••• ••• ••• 

0 0 0 ... -a 

0 0 0 ... a 

Figure 7 shows a central composite design for k=3. The value of a and the number of 
design points at the center of the design are varied to meet certain conditions. In the 
following, those conditions are chosen assuming that the approximating polynomial function 

is given by Equation (56). 

3.4. 1.1 Single center point rotatable second order experimental designs: 

A design is said to be rotatable when the variance of the estimated response--that is, the 
variance of % which in general is a function of position in the design space, is instead only 
a function of the distance from the center of the design and not on the direction. In other 
words, a rotatable design is one for which the quality of the estimator ? is the same for two 
points that are the same distance from the center of the design [3]. It is possible to develop 
central composite designs which have a single center point The value of a which will yield 
these rotatable second order designs are given in Table 3.2. 
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Table 3.2. Value of a for single center point rotatable central composite designs 


k 

a 

2 

1.414 

3 

1.682 

4 

2.000 

5 

2.378 

5 (1/2 rep) 

2.000 

6 

2.828 

6 (1/2 rep) 

2.378 

7 

3.364 

7 (1/2 rep) 

2.828 

8 

4.000 

8 (1/2 rep) 

3.364 


Note in Table 3.2 that a rotatable second order experimental design can be obtained with 
a fractional factorial design augmented with additional design points as well as with a 
augmented full factorial design. 


3.4. 1.2 Multiple center point rotatable uniform precision designs: 

In general, the variance of y varies with distance from the center of the design. However, 
by varying the number of center points, N, the variance at a distance of unity from the 
center can be made approximately equal to the variance at the center of the design. Such 
designs are referred to as uniform precision designs . The uniform precision design is based 
on the philosophy that in the central region of the design space there should be uniform 
importance as far as the variance of response is concerned, as opposed to, for example, a 
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situation in which the variance is low in the center of the design but increases drastically as 


one moves away from the design center [3]. The number of center points, m, and the value 
of a can be varied so as to obtain a rotatable uniform precision designs . Table 3.3 gives 
those values. 

Table 3.3. Values of m and a for multiple center point rotatable uniform precision designs 


k 

m 

a 

2 

5 

1.414 

3 

6 

1.682 

4 

7 

2.000 

5 

10 

2.378 

5 (1/2 rep) 

6 

2.000 

6 

15 

2.828 

6 (1/2 rep) 

9 

2.378 

7 (1/2 rep) 

14 

2.828 1 

8 (1/2 rep) 

20 

3.364 I 


3.4. 1.3 Single center point orthogonal central composite designs: 

An orthogonal central composite design can be developed where [Z]‘[Z] is diagonal. To 
obtain a design of this type a single center point can be used and the a value are taken from 
Table 3.4. 
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Table 3.4. Values of a for single center point orthogonal central composite designs 


k 

a 

2 

1.000 

3 

1.216 

4 

1.414 

5 

1.596 

6 

1.761 

7 

1.910 

8 

2.045 


3.4.1.4 Rotatable orthogonal designs: 

By varying the number of designs at the design center, m, and by selecting appropriate 
values for a, an orthogonal rotatable central composite design can be obtained. Values of 
m and a for such a design are given in Table 3.5. 
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Table 3.5. The value of m and a for multiple center point orthogonal rotatable central 
composite designs 


k 

m 

a 1 

2 

8 

1.414 1 

3 

9 

1.682 

4 

12 

2.000 

5 

17 

2.378 1 

5 (1/2 rep) 

10 

2.000 

6 

24 

2.828 

6 (1/2 rep) 

15 

2.378 

7 (1/2 rep) 

22 

2.828 I 

8 (1/2 rep) 

33 

3.364 | 


3 4 2 Discussion of the central co mposite design 

Orthogonal central composite designs have been shown to give a variance of response 
comparable to that obtained with a full 3 k factorial design. Thus, their use is justified when 
one has experimental error in the response function. Rotatable and uniform precision 
designs attempt to control the response variance. Thus there use is also justified when one 
has experimental error in the response function. However, when building a response surface 
for a deterministic model where there is no experimental error in the response function, 
their use is justified only by how well they perform of trial problems. Likewise, the designs 
were developed for the approximating function of Equation (56). If a different second order 
polynomial approximating function such as in Equation (8) were used or if a neural net was 
used to develop the response surface, then again the justification for the use of the various 
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central composite designs would have to be based on their performance on trial problems. 
Performance of various central composite designs on trial problems is next reported. 


3.4.3 Example -- Fox’s Banana Function 
Fox investigated in Reference [16] a function 

y~ 10jc 1 4 -20jc 2 x 1 2 + lQx^ +x 1 2 -2x, +5 (60) 

which has banana shaped contours as seen in Figure 8. The region of interest to be 
considered is (-1.5<x 1 <1.5, -.5<x 2 <2.0). 

A second order polynomial approximation is to be made of this function using an orthogonal 
polynomial approximation as in Equation (56). A two variable orthogonal polynomial 
approximation is of the form 

y = b g +b l x 1 +b 2 x 2 +b il ( x l - x ^ + b- n { x 2~ x ^)+b l ^c ) x 1 ( 61 ) 


where 



N 



«-i 


N 


(62) 


and where 
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( 63 ) 


N=the number of design points and 
Xj =Xj at the design points 


In the first phase of this example, Fox’s function was approximated using the second order 
orthogonal polynomial of Equation (61). The designs used in making the approximation 
were 

1. a full 5 2 factorial design, 

2. a full 3 2 factorial design, 

3. single center point rotatable central composite design, 

4. multiple center point rotatable uniform precision central composite design, 

5. single center point orthogonal central composite design, 

6. multiple center point rotatable orthogonal central composite design, 

7. minimum point design from program DESIGNS, 

8-10. minimum point design from program DESIGNS augmented by additional randomly 

selected design points, and 

11. nine randomly selected design points. 

Once an approximation was obtained, the approximate function was evaluated at a 3 1 x 3 1 
grid of points over the region of interest. The approximate function values at these 961 
points were used to develop the error parameter v G from Equation (6). Because there are 
a differing number of functional evaluations required for each of the sundry designs tested, 
a comparison of the designs based on v G is misleading. For example, the full 5 2 factorial 
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design has 25 design points each requiring a functional evaluation where as the multiple 
center point rotatable orthogonal central composite design has but 16 design points 
requiring 9 functional evaluations (in the following it is assumed that the function being 
approximated has no experimental or numerical error and thus the 8 design points at the 
design center require but one functional evaluation). Thus a comparison of performance 
based only on quality of fit is not a fair comparison. The 5 2 factorial might do a better job 
of approximating a function but the computational cost of the 25-9 = 16 extra functional 
evaluations might make it a less desirable design. 

For each design, design j, a measure of efficiency, Ej, was developed where 

£ - design j "^design ]_ (64) 

design 1 '^'design 1 

where T is the number of functional evaluations required for a given design. 

The efficiency of all the designs was compared to design 1, the 5 2 factorial design. Table 
3.6 gives, for each design tested, the number of design points, N; for central composite 
designs, the number of design points at the center of the design, m; the number of 
functional evaluations required, T; the value of v; the value of v G ; and the value of E:. 
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Table 3.6. Performance of various designs on Fox’s Banana Function, orthogonal 
polynomial approximating function, -1.5<X|<1.5, -.5<x 2 <2.0 


Design 


N 


m 


E s 


5 2 factorial design 


25 


25 


70.76 


3 2 factorial design 


64.07 


78.92 

102.46 


single center point rotatable 
central composite design 


54.36 


77.34 


multiple center point rotatable 
uniform precision central 
composite design 


13 


53.08 


77.34 


single center point orthogonal 
central composite design 


64.07 


102.46 


1.00 

.47 


.35 


.35 


.47 


multiple center point rotatable 
orthogonal central composite 
design 


16 


51.62 


minimum point design from 
program DESIGNS 


0 


77.34 


162.62 


.35 


.49 


minimum point design from 
program DESIGNS augmented 
by 2 randomly selected design 
points 


8 


8 


43.27 


105.16 


minimum point design from 
program DESIGNS augmented 
by 3 randomly selected design 
points 


53.53 


88.63 


minimum point design from 
program DESIGNS augmented 
by 4 randomly selected design 
points 


10 


10 


53.05 


86.44 


.43 


.40 


.44 


random--9 points 


21.05 


460.96 


2.10 
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Several items can be noted in Table 3.6: 


1. The design composed of 9 randomly selected design points did poorly. Even though the 
design points were chosen randomly, it turned out that the design points were not well 
scattered in the design space but were heavily concentrated in one quadrant of the design 
space. The polynomial approximation fitted the function well at the design points but poorly 
over the region of interest. 

2. The value of v G was approximately the same for the single center point rotatable central 
composite design, the multiple center point rotatable uniform precision central composite 
design, and the multiple center point rotatable orthogonal central composite design. These 
three designs differ only in the number of design points at the center of the design space. 
These designs have 1, 5, and 8 designs at the center, respectively. The effect of putting 
more designs at the center is to translate the response surface toward the center response. 
For this problem, however, the actual and approximated response were very close at the 
design center point, even for only 1 design point at the center. Thus, adding more design 
points at the design center did little to translate the response surface and thus did not 
material effect the value of v G . 

3. The eleven designs of Table 3.5 were next used to build an approximation using the 
standard second order polynomial approximation of Equation (8) instead of the orthogonal 
polynomial approximation of Equation (61). Results identical to those of Table 3.5 were 
found. The type of approximating polynomial may effect variances but does not affect 
quality of fit at the design points or over the region of interest. For those problems were 
there is no experimental or numerical error associated with functional evaluations, one is 
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not interested in variance. Thus, there is little advantage in using the orthogonal polynomial 
approximating functions over a standard second order polynomial function. 

4. Based on efficiency, the single center point rotatable central composite design, the 
rotatable uniform precision central composite design, and the rotatable orthogonal central 
composite design performed the best but none of the designs gave a good approximation 
over the region of interest. Over a small region of interest, one could expect that a second 
order polynomial approximation could well approximate the given function. Obviously, here 
the region of interest is too large for a second order approximation to be a good one. Thus 
a smaller region of interest was chosen, -.5 <x lv 5, -.5 <x 2 < .5. Table 3.7 compares the eleven 
designs using this region of interest. Notice that over this smaller region of interest, all the 
designs gave a much better approximation to the function. 

5. For the smaller region of interest, based on efficiency, the 3 2 factorial design, the single 
center point orthogonal central composite design, and the augmented minimum point 
designs performed the best. Obviously, the optimum choice of design is problem dependent. 
However, all designs except the randomly selected design performed much better than the 
5 2 factorial design. 
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Table 3.7. Performance of various designs on Fox’s Banana Function, orthogonal 
polynomial approximating function, -.5<X!<.5, -.5<x 2 <.5 


Design 

N 

m 

D 

— 

V G 

M 

5 2 factorial design 

25 

m 

25 

11.16 

8.57 

1.00 

3 2 factorial design 

9 

... 

9 

13.27 

10.95 

.46 

single center point rotatable 
central composite design 

9 

i 

9 

6.58 

14.74 

.62 

multiple center point rotatable 
uniform precision central 
composite design 

13 

5 

9 

5.88 

14.74 

.62 

single center point orthogonal 
central composite design 

9 

1 

9 

13.27 

10.95 

.46 

multiple center point rotatable 
orthogonal central composite 
design 

16 

8 

9 

5.47 

14.74 

.62 

minimum point design from 
program DESIGNS 

6 

■ 

6 

0 

18.66 

.52 

minimum point design from 
program DESIGNS augmented 
by 2 randomly selected design 
points 

8 

... 

8 

5.74 

11.82 

.44 

minimum point design from 
program DESIGNS augmented 
by 3 randomly selected design 
points 

9 

... 

9 

6.45 

10.53 

.44 I 

minimum point design from 
program DESIGNS augmented 
by 4 randomly selected design 
points 

10 

... 

10 

6.33 

10.29 

.48 

random-9 points 


■ 

9 

2.42 

47.22 

1.98 
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3.4.4 Conclusion 

Second order polynomial approximations or neural net equivalents are often adequate for 
building response surfaces, especially if the region of interest is small. Central composite 
designs are convenient for building the second order approximations. They provide the 
necessary information for determining all of the coefficients of the approximating polynomial 
and give a good distribution of points in the design space. The approximating function can 
be made to closely fit the exact function at the design center by using multiple center points. 
When modeling deterministic systems, each functional evaluation at the design center yields 
the same function value. Thus, for deterministic models, only one functional evaluation 
need be performed at the center point even when multiple center points are used. Table 
3.8 gives information relevant to central composite designs for various number of design 
variables, k. Central composite designs give over-determined second order polynomial 
approximations. In other words, there are more design points in the design than there are 
undetermined coefficients in a second order polynomial approximation. Table 3.8 also gives 
the percentage that the approximation is over-determined. Previous studies [19] have 
indicated that designs which give approximations that are around 20-50% over-determined 
tend to be efficient designs. One can see that the central composite designs are reasonable 
for k<6. For larger k values, too many design points are being used by the central 
composite designs. For k>5, an augmented minimum point design is a better choice. 
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Table 3.8. Information relevant to central composite designs for various number of design 
variables 


Number of design 
variables, k 

Number of 
coefficients in a 
2nd order 
polynomial 
approximation 

Number of 
functional 
evaluations 
required with a 
central composite 
design 

% over-determined 

1 

3 

4 

33 

2 

6 

8 

33 

3 

10 

14 

40 

4 

15 

24 

60 

5 

21 

42 

50 

6 

28 

76 

171 

7 

36 

142 

294 

8 

45 

272 

504 
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4. Optimality Criteria 


4.1 D. A. E. G. and V O ptimality Criteria 

It was pointed out in Section 3 that even for a small number of factors, a complete factorial 
experiment become impractical if functional evaluations are computationally or 
experimentally expensive to obtain and thus one is forced to use some sub-set of the 
factorial design or an alternate design requiring fewer experiments. Section 3 shows that 
the variances of the coefficients of a polynomial approximation and the variance of the 
predicted response involve the matrix [H] given in Equation (33) and repeated here: 

m-craTED" 1 (65) 

Schoofs [22] lists five criteria for selecting a sub-set of the factorial designs. These criteria 
involve the matrix [H]. The criteria, referred to as optimality criteria, attempt to make [H] 
minimal. However, "the minimum of a matrix is not a well defined concept and a number 
of operational criteria have been developed" [22]. The optimality criteria for selecting a 
subset of a full factorial design can be based on selecting the subset satisfying the following 
criteria: 

1. D-optimality, which is achieved if the determinant of [H] is minimal which in term gives 
that the product of the eigenvalues of [H] is minimal. 

2. A-optimality, which is achieved if the trace of [H] is minimal which in term gives that 
the sum of the eigenvalues of [H] is minimal. 

3. E-optimality, which is achieved if the largest eigenvalue of [H] is minimal. 
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4. G-optimality, which is achieved if the maximum over all candidate points of the 
estimated response variance is minimal. 

5. V-optimality, which is achieved if the estimated response variance, averaged over all 
candidate points is minimal. 

4.1.1 Criteria Applied to a One Dimensional Example 

An example is considered here to compare the performance of the 5 optimality criteria. 
The following test function of one variable was considered: 

y=2+x+sin[— (x+ 1)], -lsxsl (66) 

2 


This function was approximated with polynomials of order 1-4. The approximations shown 
in Figure 9 were developed using 13 designs, uniformly spaced in the region of interest. 
These approximations were then used to generate the functional values at 61 uniformly 
spaced points in the region of interest which were used to plot the curves of Figure 9. 

Further approximations of Equation (66) were developed using various number of design 
points, n. The designs selected were 

1. uniformly spaced design points, n= 5,7,9, 11, 13; 

2. randomly selected design points, n= 5,7,8, 11, 13; 

3. an n member subset of the 13 uniformly spaced design points, n= 5,7,9, 11. 
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Under item 3, the subset of design points was chosen using: 

1. D-optimality, 

2. A-optimality, 

3. E-optimality, 

4. G-optimality, and 

5. V-optimality. 

A FORTRAN program was written to perform the investigation under item 3. The 
demanding part of the programming was to identify all the possible subsets from the set of 
thirteen design points. After developing a procedure to identify all combinations, each 
subset was used to build the [H] matrix. The "optimal" [H] matrix was then determined 
using the five optimality criteria. The coefficient v G was then computed for the optimal 
subset. Figures 10-13 show the value of v G for the D, A, E, and G optimality criteria when 
a first, second, third, and fourth order approximation is being made, respectively, versus the 
number of design points specified in the subset . Also shown in those figures is the value 
of v G for designs consisting of design points uniformly spaced in the region of interest. 

It was found that for all subsets of size r from a design point set of size n that the estimated 
response variance, averaged over all candidate points, was invariant. This finding 
undoubtedly could also be proven theoretically but such a proof was not attempted. From 
this example, one can conclude that the V optimality criteria, which employees the estimated 
average response variance, is not a viable criteria for selecting a subset of design points from 
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a given set- From Figures 10-13, one can see that in most cases there is little difference in 
the performance of the various optimality criteria with criteria D and G performing slightly 
better than the other two criteria. As can be seen in Figure 12, on one occasion (when 
using a third order polynomial approximation and when selecting a subset of 5 design points 
from the 13 design point set) the G optimality criteria performed poorly while the D criteria 
did not. Thus, this example indicates that the D optimality criteria mav be the criteria of 
£hoi££. There is a further advantage in using the D optimality criteria. The requirement 
that the determinant of [H] is minimal is equivalent to a requirement that the determinant 
of [G] is maximal where 

[G]=[Z]1Z} (67) 

Thus the D optimality criteria insures that the procedure for determining polynomial 
coefficients in Equation (12) will be well defined. In other words. Equation (12) uses the 
inverse of [G]. The D optimality criteria guarantees that [G] is not singular. 

One can see in Figures 10-13 that, in most cases, all the optimality criteria performed worst 
than the uniformly spaced design case. This example indicates that a design picked using 
an optimality criteria m ay be no bet ter than a design of the same size in which the design 
points are uniformly located in the design space 

4.2 S and O Optimality Criteria 

The previous optimality criteria involved only the matrix [H] and did not consider the 
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function to be approximated. Thus for a given number of design variables and level of 
approximation, the same designs would be selected no matter what the nature of the 
function to be approximated. Initially it was thought that a superior optimality criteria 
would have to consider the nature of the function. Thus two additional optimality criteria 
were examined: 

1. S-optimality, which is achieved if the average error of approximation at the design points 
is minimal and 

2. Q-optimality, which is achieved if the maximum error of approximation at the design 
points is minimal. 

Here 

^(yr9i> 2 ( 68 ) 

average error of approximation^ 


and 


maximum error of approximation=mzx (y^y) 2 , i- h-jr 


(69) 


where r is the size of the subset of design points to be selected. One can see that with the 
S and Q optimality criteria, the function to be approximated effects the design points 
selected. 
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4,2.2 Criteria Applied to a One Dimensional Example 

The one dimensional example problem of Section 4.1.1 was then re-examined. Figures 14- 
17 show values of v G using the S and Q optimality criteria and using a first, second, third, 
and fourth order polynomial approximation, respectively, versus size of the subset of design 
points. Also shown in these figures are results for uniformly spaced design points. One can 
see in these figures that terrible approximations were obtained with these criteria when only 
small subsets of design points were selected from the original set. Figures 18-20 indicate 
why such bad approximations are obtained with these two criteria. 

Figure 18 depicts results obtained by having eleven designs points selected, using the Q 
optimality criteria, from a set of 13 design points. The Q optimality criteria finds an 
approximation such that the maximum error of the approximation over eleven design points 
is minimal. One can see in Figure 18 that the approximating function did indeed well fit 
the exact function at the 1 1 design points selected. However, the approximating function 
did a poor job of approximation at the ends of the region of interest and thus would not 
yield a low value of v G . Figure 19 is similar to Figure 18 except that this figure depicts 
results obtained by having 7 design points selected from the set of 13 design points. One 
can see that for the optimum design selected, there is an almost perfect approximation at 
the design points selected but over a much larger region the approximation is poor and thus 
a large value of v G would be obtained. In Figure 20, only 5 design points are selected. 
Again at those design points, an almost perfect approximation is obtained but a terrible 
approximation is obtained over a large part of the region of interest and thus a large v G 
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would be obtained. Thus we can conclude that the S and Q optimality criteria are not 
operative . 


4.3 An Alternate Approach-Rand om Selection of Designs 

The effect of randomly picking design points was next considered. Here designs are picked 
in the region of interest using a random number generator. 

4 31 Random Selection of Designs Applied to a One Dimensi onal Example 
For the one dimensional problem under consideration, first, second, third, and fourth order 
approximations were considered. Design point sets containing 5,7,9,11, and 13 design points 
were developed by randomly picking design points in the region of interest using a random 
number generator. Approximations were developed using the design sets. Results using 
these approximations are compared in Figures 21-24 to results using uniformly spaced design 
points. One can see in these figures that most of the time results from randomly picked 
design points are either as good as or not much worst than results from uniformly spaced 
design points. However, on two occasions, when the number of design points in the design 
set was small, a relatively poor approximation was obtained. Obviously where one is picking 
only a small number of points using a random number generator, there is a chance that a 
bad set of points can be generated and indeed on these two occasion a poor selection of 
points was made. In general however, when more design points are randomly selected, 
those points should be scattered throughout the design space and good approximations 
should be obtained, in conclusion, randomly select ing desi g n p oints may be a viable metho d 
of design selection . 
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4.4 Larger Problems 

Consider a problem in two variables and consider that the potential design points will be 
taken from a 6 x 6 grid of points. Let 


r = total number of design points in the set of potential design points, 
c = number of design points in the selected subset of design points, 
nc = the number of different combinations of designs in the subset. 


For the problem at hand, r =36. Subset sizes of c = 15, 20, 25, and 30 are to be considered. 
The number of possible combinations of design points in the subset, nc, is given by 


nc= 


rl 

(r-c)! c! 


(70) 


Table 4.1 summarizes the number of combinations for this study. 


Table 4.1 Number of combinations of designs in a two variable study 


r 

Total number of design 
points 

c 

Number of point in subset 

nc 

Number of combinations 

36 

15 

5,567,902,560 

36 

20 

7,307,872,110 

36 

25 

600,805,296 

36 

30 

1,947,792 
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One can see that for even small problems, it is infeasible to examine all possible 
combinations of subsets of size N from a given set of design points. Welch [23], instead of 
evaluating all possible N-point designs, developed a "branch and bound" algorithm which 
guarantees global D-optimal designs but which does not generate and evaluate all possible 
designs. However, even here the computing costs are high. Fedorov [24] developed another 
technique which neglects the integer character of the components of the design set and 
obtains a discrete design which is rounded off to an exact design. Reference [22] reports 
that these designs are considered only approximate. The most popular algorithm seems to 
be DETMAX by Mitchell [25]. Quoting reference [22], "The algorithm starts with an initial 
m-point ED (experimental design); the final goal is an optimal N-point ED. During each 
iteration step that candidate point, which results in the largest increase of det(M), is added 
to the design, and subsequently that point, which results in the smallest decrease of det(M), 
is removed from the design. The number m of points in the initial design may be larger or 
smaller than N. If necessary the algorithm first adds (if m<N) or rejects (if m>N) points 
until the number of points in the ED is equal to N. In order to avoid local optima the 
algorithm is able to perform ’excursions’, in which several points are added at one go and 
subsequently the number of points is reduced to N. If the resulting N-point ED has not 
been improved, another excursion will be made from the same initial design. If the 
excursion is successful the resulting ED will be used as starting ED in a further attempt to 
maximize det(M). The algorithm terminates when, after several excursions, no better ED 
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is found. The algorithm generates high quality EDs against relatively low computing costs." 
An attempt is being made to obtain the algorithm DETMAX. 

4.5 Op timality Criteria Based on Minimizing Uncertainty 

Reference [21] considers problems where there is no experimental error. That reference 
uses an optimality criteria based on selecting a design which minimizes the uncertainty in 
the approximating function. That reference was given mixed reviews by a number of leading 
authorities in the field [21] (reviews follow the paper). The formulation is quite theoretical 
and difficult to follow. The formulation seems to have promise but requires additional 
theoretical development before it becomes operative. 

4.6 Conclusion 

There is little rational for using any of the investigated optimality criteria when building 
approximations of functions which contain no experimental error. However, the D- 
optimality criteria can conveniently be used as a heuristic in selecting design points. 


Previous investigations have indicated that approximations should be over-determined. That 
is to say that more training pairs should be used to build an approximations than the 
number of associated undetermined parameters. It has been suggested that a 20-50% over- 
determined system might be reasonable. The program DESIGNS, described in Section 2, 
develops enough designs to exactly determine a quadratic approximation of a given function. 
The D-optimality criteria can be used as a heuristic for selecting design points to 
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supplement those generated by DESIGNS. The use of the D-optimality criteria to select 
the supplementary points would guarantee than no singular matrices would be encountered 
in determining the undetermined parameters associated with the polynomial approximation. 


56 


5. Significance Testing of Coefficients 


5.1 Introduction 

When the training pairs used to build a polynomial response surface contain experimental 
or numerical error, certain coefficients in the polynomial approximation may not be 
significant. In other words, even though one calculates a value for some coefficient, b it the 
experimental or numerical error may have such an effect on that coefficient that it could just 
as well be taken as zero as the value calculated. In situations like this, it may be 
advantageous to drop that term from the polynomial approximation and redevelop the 
response surface. Such a procedure is discussed in pages 34-38 of [3] and an automated 
procedure for performing such an operation was developed in [26]. Testing of significance 
involves the t-test which is next described. 

5.2 t-test 

Coefficients of the polynomial approximation are found from Equation (12). The 
determination of those coefficients involve the matrix [H] where 

m-mmr' (71) 


A number of terms must now be defined: 


N 


E 


mean square error =MSE- 


i-i 


N-m 


(72) 
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standard error coefficient=se i = s jMSE H u 


(73) 


'r 



(74) 


where 

N = the number of design points and 

m = the number of coefficients in the polynomial approximation. 

In making the test of significance, t ; from Equation (74) is compared to tabulated values of 

t a . The value of t a is taken from a table of "Percentage Points of the Student’s t 

Distribution" [3). The value taken depends on the level of significance desired. In lieu of 

using tabulated values, t a is often taken as four [26]. If t, is less than t a (tj < t a ), then that 

0 

coefficient’s importance in approximating the response is deemed to be insignificant and 
therefore may be eliminated from the response function. 

The primary focus of this study was to examine methods of developing good response 
surfaces for deterministic models, i.e. for systems that contain no experimental or numerical 
error. Statistical testing of coefficients presupposes experimental or numerical error and 
thus is not relevant when approximating response which contains no error. However, the 
method was thought to perhaps offer a heuristic for improving the quality of a response 
surface even if experimental or numerical errors are not present. Thus, two examples were 
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examined. Results are next reported. 


5.3 Example 1 - Fox’s Banana Function 

Example 1 again examines Fox’s Banana Function [16]. A complete second order 
polynomial approximation (m = 6) and a complete third order polynomial approximation 
(m=10) were developed. These approximations were developed using a complete 6 2 
factorial design (N = 36). A t-value, t,, was calculated for each parameter, b i; and compared 
to t a = 4. Parameter that lack significance (t, < t a ) were eliminated. A new approximation 
was then developed using only the significant parameters. The values of v and v G from 
Equations (5) and (6), respectively , were developed for the complete polynomial and for 
the polynomial containing only terms deemed significant. Results are shown in Figures 25 
and 26. On can see in these figures that eliminating coefficients deemed insignificant had 
an adverse effect on the quality of the approximation over the region of interest. 

5.4 Example 2 

The effect of eliminating coefficients deemed insignificant was tested on the function 

T=(4+x 1 ) 3 +sin[^(x 1 +l)]+2+x 2 4 +sin(-)+7x 2 x I (75) 

2 2 


Again, a complete second order polynomial approximation (m=6) and a complete third 
order polynomial approximation (m = 10) were developed. These approximations were 
developed using a complete 6 2 factorial design (N=36). A t-value, tj, was calculated for 
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each parameter, bj, and compared to t a = 4. Parameter that lack significance (t ; < t a ) were 
eliminated. A new approximation was then developed using only the significant parameters. 
The values of v and v G from Equations (5) and (6), respectively , were developed for the 
complete polynomial and for the polynomial containing only terms deemed significant. 
Results are shown in Figures 27 and 28. On can see in these figures that eliminating 
coefficients deemed insignificant offered no improvement in the quality of the response 
surface. 

5.5 Conclusion 

The applicability of significance testing of polynomial coefficients when modeling 
deterministic systems was considered. Two examples were examined to see if eliminating 
terms of polynomial approximations which were deemed to be insignificant by the t-test 
would improve the quality of the response surfaces developed. Based on these two 
examples, it was concluded that no improvement in the predictive capability of response 
surfaces over regions of interest would be obtained with such a procedure. The relevance 
of significance testing is when modeling systems containing numerical or experimental error. 
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6. Applicability of the Response Surface Technique 


6.1 Introduction 

The following study was performed to ascertain under what circumstances could the 
response surface technique be used to advantage in engineering optimization application. 
In this regard, assume that a quadratic polynomial approximations is to be made of functions 
of n variables. The number of undetermined coefficients in that approximation is: 


number of coefficients = 


(n+l)(n+2) 

2 


(76) 


Previous studies [19] have shown that the best approximations are obtained when the 
approximations are over-determined. Thus, the number of functional evaluations required 
to make the approximation is: 

number of functional evaluations- - ^ n+ ( 77 ) 


where S determines the degree that the approximation is over-determined. 

The functional evaluations required to build the approximation are initially performed 
before the start of the optimization process. By using parallel processing, these functional 
evaluations may be less computationally expensive than evaluations made sequentially in a 
direct optimization procedure. The number of required evaluations of Equation (77) is then 
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equivalent to a reduced number of sequential evaluations thus. 


equalivent number functional evaluations = 


6P(n+l)(n+2) 

2 


(78) 


where /3 is a coefficient of efficiency associated with parallel processing. 

An optimum solution can be attempted using the response surfaces developed instead of the 
original functions. However, because of the inexact nature of the approximations, a new set 
of response surfaces may have to be developed at the most recent approximate solution and 
another optimal solution attempted. This procedure may have to be repeated a times to 
reach the optimum solution for the original problem. The total number of equivalent 
f un ctional evaluations performed in reaching this optimum is: 

total equivalent functional evaluations = — (79) 


If the solutions was attempted by direct optimization techniques instead of using response 
surfaces, Barthelemy [27] states that a solution can be obtained in most cases using no more 
than ifr first derivative evaluations. If the first derivatives are obtained by finite difference 
formulae, an estimate of the number of functional evaluations required by a direct solution 

procedure is: 

functional evaluations direct methods*1f(n+ 1) 
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If the response surface technique is to be competitive with the direct solution technique, 
then from Equations (4) and (5) one must have: 


gp6(n+l)(n+2) 

2 


£yi|r(»+l) 


(81) 


where y is a convenience factor associated with using response surfaces. In other words, an 
investigator may tolerate more functional evaluations with the response surface technique 
than with the direct solution procedure just for the convenience of using response surfaces. 
Rearranging Equation (81) gives 

[n+l][ gp6(/,42) -»]*0 (82) 

2y 


Since (n+1) is positive, one obtains 

ttP6(n+2) 

2y 


(83) 


or 


n^-2 
a p6 


(84) 
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In review 


a = number sequential optimizations 
P= parallel processing coefficient 
6 = overdetermined coefficient 

y = convenience coefficient 
i|r= direct solution coefficient 


Reasonable ranges of the parameters involved are 


a =1.00-4.00 
p =0.10-1.00 
5 = 1.25-1.75 
Y = 1.00-3.00 

ifr =6.00- 10.0. 


( 86 ) 


For an approximate upper bound on the number of design variable that could be 
economical used with the response surface technique take: 


a = 1.00 
p=0.10 
6 = 1.25 
y =3.00 
\|r = 10.0 


(87) 


giving 


n<;498 


( 88 ) 


Under the most unfavorable set of circumstances, that is: 
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( 89 ) 


a =4.00 
P = 1.00 
5 = 1.75 
Y = 1.00 
i|r=6.00 


one obtains 


n-0 


(90) 


Thus depending upon the problem, one could use the response surface technique for n = 0 
to n=500 variables. Consider the following reasonable set of parameters 


<x=3.00 
P =0.50 

6 = 1.25 (91) 

y = 1.50 
= 8.00 


giving 


13 (92) 


Thus, it is reasonable to assume that the response surface technique could be used for up 
to 10-15 design variables. 
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6.2 Conclusion 


Under the most advantageous circumstances, the response surface technique applied to 
engineering optimization application could be used for up to 500 design variables. Under 
the worst set of circumstances, it is entirely inappropriate. Under normally expected 
circumstances, this technique might be used to advantage for 10-15 design variables. 
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7. Additional Examples 


7.1 Introduction 

The next several examples examine the effect of design selection on the quality of 
approximations. In each case, a second order polynomial approximation is made of a trial 
function. Different number of design variables are considered in each example. Thus, for 
each example different designs are appropriate. In the first example, there are 4 design 
variables. When there are fewer than 6 design variables, central composite designs are a 
possible appropriate choice. Other choices are the 3 k factorial design, the minimum point 
design, the augmented minimum point design, or randomly selected design. All of these 
designs are considered in that example. In the second and third examples, there are 15 and 
20 design variables, respectively. Here, the 3 k factorial design and central composite designs 
contain too many design points to be practical. For these examples, the minimum point 
design, the augmented minimum point design, and the randomly selected design are 
appropriate and are considered. 

7.2 The 35 Bar Truss with 4 Design Variables 

In many response surface applications, the function to be approximated is a relatively 
smooth function of the design variables which can be approximated with a lower order 
polynomial or an artificial neural net with only a few nodes on the hidden layer. A problem 
of this type is shown in Figure 29. In this example, all loads shown in Figure 29 are in kips, 
all members of the lower chord of the truss are assumed to have area, A v and all members 
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of the upper chord to have area, A 2 , all vertical and diagonal members to have area, A 3 . 
The depth of the truss is H. A response surface is to be constructed for the stress in 
member BC in terms of the design variables, Xj thus 


x t -l/A p i= 1,3 
x 4 =. 09375#- .4375 


(93) 


The region of interest is 


2 in 2 8 in 2 

6 ft zHs 10ft 


(94) 


or in terms of the design variables 

.125 ^.5 


(95) 


A number of designs were used to develop a second order polynomial approximation for the 
stress in member BC. Each approximation was then used to predict stress on a 5 x 5 x 5 
x 5 grid of points. The predicted stress and the actual stress on these NG = 625 grid of 
points were then used to develop v G from Equation (6). The parameter v G is a measure of 
the quality of the approximation over the region of interest. 


The different designs examined required different numbers of functional evaluation. So as 
to get a measure of the quality of fit of the approximation over the region of interest which 
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Table 7.1 The 35 bar truss with 4 design variables, 2nd order polynomial approximation 


Description 

m 

a 

ID 


v(%) 

v G (%) 

Ei 

3 4 factorial design 

... 

... 

81 

81 

3.34 

2.41 

1.00 

single center point rotatable 
central composite design 

i 

2.000 

25 

25 

0.66 

2.67 

0.34 

multiple center point rotatable 
uniform precision central 
composite design 

■ 

2.000 

31 

25 

0.59 

2.67 

034 

single center point orthogonal 
central composite design 

i 

1.414 

25 

25 

1.47 

237 

0.30 

multiple center point rotatable 
orthogonal central composite 
1 design 

12 

i 

2.000 

36 

25 

0.55 

2.67 

0.34 

minimum point design from 
program DESIGNS 

n 

... 

15 

15 

0.00 

3.99 

0.31 

minimum point design from 
program DESIGNS augmented 
by 3 randomly selected design 
points 

it# 

... 

18 

18 

0.40 

3.86 

0.36 

minimum point design from 
program DESIGNS augmented 
by 6 randomly selected design 
points 

••• 

... 

21 

21 

0.38 

3.91 

0.42 

minimum point design from 
program DESIGNS augmented 
by 9 selected design points 

... 

... 

24 

24 

0.41 

3.77 

0.46 

randomly selected design 

... 

... 

25 

25 

0.00 

824.2 

105 


m= number of design points at the center of the design space 
T = the total number of design points 
F = the number of functional evaluations required 
a = parameter which defines location of certain design points 
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takes into account the number of functional evaluations performed, the efficiency, Ej, from 
Equation (64) was developed for each design. Table 7.1 reports for each design considered, 
the efficiency, Ej, as well as other relevant information. 

One can see in Table 7.1 that all the designs considered, except the randomly selected 
design, gave a good approximation over the region of interest. Randomly selected designs, 
which often work well, can sometimes suffer from the problem that the coefficient matrix 
used to solve for the approximation’s associated parameters is poorly conditioned or that 
the design points selected are not well scattered throughout the design space. In either case, 
they can yield a poor approximation over the region of interest as in this example. 

The 3 4 factorial design well approximated the trial function. However, because it uses so 
many design points its efficiency measure is poor and thus is not a design of choice. The 
single center point orthogonal central composite design and the minimum point design from 
program DESIGNS performed the best, based of their efficiency. However, excluding the 
randomly selected design and the 3 4 factorial design, all of the designs considered gave a low 
value of v G and had approximately the same value of efficiency. 


Under normal circumstances, information is not available to calculate v G and one must use 
the parameter v as a measure of the quality of fit over the region of interest. However, the 
parameter v is only a measure of quality of fit over the region of interest if the 
approximation is over-determined. Thus, under normal circumstances one would not want 
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to use the minimum point design. This example indicates, that for problems of the size of 
this example, that any of the central composite designs or the augmented minimum point 
designs would be appropriate. 

7.2 The 35 bar truss with 15 design variables 

This example again considers the 35 bar truss of Figure 29. In this example, H is 10 ft., the 
areas of the 14 bars of the top and bottom chords are i = 1,14, and the area of the 
vertical and diagonal members is A 15 . The design variables of the problem are taken as 

XfUAp *=U5 (96) 


The region of interest is 


2 in 2 s A i s; 8 in 2 


(97) 


or in terms of the design variables 

.125^<;.5 (98) 

Response surfaces were developed for the stress in member BC using a 2nd order 
polynomial approximation. The approximation were developed using various designs. To 
test the quality of the approximations over the region of interest, the function and the 
approximations were evaluated at NG = 500 randomly selected test points over the region 
of interest. That information was then used to calculate v G from Equation (6). The random 
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number generator used to develop design points uses, in generating its numbers, an initial 
seed parameter, IFLAG. A different value of IFLAG was used to generate the 500 test 
points than was used to generate random points in the randomly selected designs or in the 
augmented minimum point designs. Thus, the test set of points does not duplicate any of 
the design points in the designs considered. Results of this investigation are reported in 
Table 7.2. 

One will notice in Table 7.2 that only minimum point designs, augmented minimum point 
designs, and randomly selected designs are considered. A 3 15 factorial design contains over 
14 million design points. Thus, the use of the 3 15 factorial design is out of the question. For 
a problem in k design variables, the central composite design uses a 2 k factorial design 
augmented by 2k +1 additional design points. Thus, such a single center point central 
composite design for this problem contains 32,799 design points. Flere again, such a design 
is impractical. One can develop a central composite design by augmenting only a fraction 
of the 2 k factorial design. For this problem, a single center point central composite design 
using only a 1/4 fraction of the 2 15 factorial design would contain 8,223 design points which 
is still an impractical design. Thus, for problems of the size of this example, only the 
minimum point designs, augmented minimum point designs, and randomly selected designs 
are of reasonable size. 

We can see in Table 7.2 that all of the designs with the exception of the "randomly selected- 
-exactly determined design" did a good job of approximating truss behavior. A singular 
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matrix was encountered in Equation (10) for the randomly selected-exactly determined 
design. With completely randomly selected designs, there is always the possibility of having 
a poorly conditioned coefficient matrix [Z] in Equation (10) and indeed this occurred in this 
problem. However, there was no problem with matrix conditioning using randomly selected 
over-determined designs. 


Table 7.2 The 35 bar truss with 15 design variables, 2nd order polynomial approximation 


Description 

F 

v % 

v G % 

Ei 

minimum point design from 
program DESIGN- 
exactly determined 

136 

0 

1.263 

1.0 

augmented minimum point 
design— 20% over- 
determined 

163 

0.083 

0.294 

0.28 

augmented minimum point 
design— 40% over- 
determined 

190 

0.087 

0.060 

0.07 

randonj selection-exactly 
determined 

136 

* 

* 

* 

random selection-20% 
over-determined 

163 

0.003 

0.029 

0.03 

random selection-40% 
over-determined 

190 

0.003 

0.010 

0.01 


* singular coefficient matrix 


The efficiency parameter, E jt is calculated in Table 7.2 but it is rather a meaningless 
parameter for this problem because all the designs so well fit the exact function. In real life 
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situations, one usually does not have available information for calculating v G . Thus, the 
parameter v or like term must be used as a measure of the quality of the approximation. 
The parameter v is not a meaningful measure of the quality of fit over a region of interest 
unless the system is over-determined. Thus for this example, the design of choice would be 
either the 20% over-determined minimum point design or the 20% over-determined 
randomly selected design. 


wnntrai 


[TTifnijiT 


This example considers a problem with even more design variables. The function tested is: 


20 20 20 20 20 

y= i- + E *< + E E wE E 

i«l i» 1 j-i i-1 j m i 


(99) 


A second order polynomial function was used to build the response surface approximating 
this function. The polynomial approximating function had 231 undetermined coefficients. 
Because of the large size of this problem, factorial designs and central composite designs 
are not appropriate. A minimum point design, augmented minimum point designs, and 
randomly selected designs were considered. Values of the test function and approximate 
function were evaluated at NG = 1000 randomly selected points and the parameter v G was 
developed using this information. The measure of efficiency of the designs examined along 
with other relevant information is given in Table 7.3. 
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Table 7.3 Analytical function with 20 design variables, 2nd order polynomial approximation 


Description 

F 

v % 

v G % 

E. 

minimum point design from 
program DESIGN- 
exactly determined 

231 

0 

88.93 

1.0 

augmented minimum point 
design-20% over- 
determined 

277 

5.83 

49.82 

0.67 

augmented minimum point 
design— 40% over- 
determined 

323 

9.58 

18.03 

0.28 

random selection— exactly 
determined 

231 

* 

« 

* 

random selection-20% 
over-determined 

277 

0.61 

7.21 

0.10 

random selection-40% 
over-determined 

323 

0.46 

1.20 

0.02 


* poorly conditioned coefficient matrix 


Just as in Example 7.2, a exactly determined randomly selected design gave a poorly 
conditioned coefficient matrix. These examples indicate that randomly selected exactly 
determined designs should be avoided. The 40% over-determined randomly selected design 
did an excellent job of modeling the test function and was the most efficient design 
considered. It seems that on problems with a large number of design variables that 
randomly selected over-determined designs should be expected to work well. 
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7.4 Conclusion 

The examples of this section have shown that design selection depends on the number of 
design variables. If the number of design variables is less than 6, appropriate designs are. 

1. augmented minimum point designs 

2. central composite designs 

3. over-determined randomly selected designs. 

When there are more than 6 design variables, the central composite designs contain too 
many design point for consideration. For more than 6 design variables, appropriate designs 
are then 

1. augmented minimum point designs 

2. over-determined randomly selected designs. 

The example examined indicate that in all cases, over-determined designs should be used. 
They the most efficient designs. Also, when a design is over-determined the coefficient v 
can be used as a measure of the quality of the approximation over a region of interest. 
Being able to use v as a measure of the quality of fit over the region of interest is very 
important because, in general, information is not available to determined the parameter v G . 
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8. Augmented Minimum Point Designs 


8.1 Introduction 

Design selection in the literature concentrates of linear or quadratic response surfaces. This 
study has also concentrated on quadratic approximations for several reasons: 

1. linear approximations, in most instances, will be inadequate to model functions of 
interest, 

2. for many problems, a 2nd order approximation will be adequate to model response 
especially if the region of interest is limited, 

3. there is a scarcity of literature which address design selection for cubic or higher order 
polynomial approximations, and 

4. in optimization process using response surfaces, for moderate size problems, it is more 
computationally efficient to perform a sequence of quadratic approximations than one cubic 
or higher order approximation. This fact is next discussed. 

The number of terms in a second order polynomial in n design variables is 

number terms quadratic =(n+l)+?^-Q- (100) 


The number of terms in a 3rd order polynomial in n design variables is 
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(101) 


3 ft! 

number terms cubic = 1 +-/»(/»+ 1)+— — — 

2 6(n-3)! 

Table 8.1 gives, for various number of design variables, the number of terms in a 2nd order 
and 3rd order polynomial and their ratio. 


Table 8.1 Number of terms in a 2nd and 3rd order polynomial and their ratio 


number of design 
variables, n 

number of terms 
in quadratic 

number of terms 
in cubic 

cubic/quadratic 

3 

10 

20 

2 

6 

28 

84 

3 

9 

55 

220 

4 

12 

91 

455 

5 

15 

136 

816 

6 


One can see that for problems with more than 6 design variables, it will probably be more 
computationally efficient in an optimization algorithm to utilize a sequence of quadratic 
response surfaces than one 3rd or higher order response surface. When there are 6 or fewer 
design variables, 3rd or 4th order response surfaces may be used to advantage. 

In this report, the term "minimum point design" refers to a design that has just enough 
design points to allow the determination of coefficients of an approximating polynomial. 
The term "augmented minimum point design" is a minimum point design which contains 
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additional design points. Thus, augmented minimum point designs are over-determined 
designs. The studies that have been performed in this report indicate that augmented 
minimum point designs are competitive with, if not better than, central composite designs 
for developing a 2nd order response surface. A program DESIGNS [20] was developed for 
generating augmented minimum point designs for developing a 2nd order response surface. 
That program is described in Section 8.2. 

When there are 6 or fewer design variables, it may be computationally beneficial to use a 
3rd order or 4th order response surface. Thus, the program DESIGN4 [28] was developed 
to generate augmented minimum point designs for a 4th order response surface. The 
program DESIGN4 is discussed in Section 8.3. The program can also be used to develop 
a 3rd order response surface. The 3rd order minimum point design is a subset of the 4th 
order minimum point design. Thus the 4th order minimum point design will give an over- 
determined 3rd order approximation. Additional randomly selected design points can be 
added to the 4th order minimum point design to give the desired degree that the 3rd order 
approximation is to be over-determined. 

8.2 Augmented Minimum Point Designs for 2nd Order Approximatinnc 


The basic building block for program DESIGNS is the star pattern of design points. Figure 
4 shows the star pattern for 3 design variables. This pattern of design points allows one to 
determine those coefficients of a 2nd order polynomial approximation associated with the 
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terms 


1, x t , xf, i= Ijt 


( 102 ) 


To be able to determine the coefficients associated with the terms 

x ?r '*i 


(103) 


one must supplement the star pattern with one additional design point in the x^ planes. 
Figure 30 shows the additional design point in the Xj, Xj plane. Figure 6 shows the total 
minimum point design for 3 design variables. 

Studies of this report indicate that designs should be over-determined. Having a design that 
is 20-50% over-determined is a good compromise between keeping down the number of 
design points while still getting a good approximation. The program DESIGNS augments 
the minimum point design with a user selected number of random design points. 

8.2.1 Specifics of program DESIGNS 

A listing of the FORTRAN program DESIGNS is found in Appendix 1 and a copy of that 
program is found in file "designs.f on the floppy disk accompanying this report. The 
program should be compiled with a F77 compiler with the compiled program called "design". 
To run the program just enter "design" from the keyboard. The program prompts the user 

for 
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1. the number of design variables, 

2. the number of designs points to augment the minimum point design, and 

3. a seed parameter, IFLAG, which is used to generate the random numbers (IFLAG can 
be entered as any positive integer). 

The program then generates a design in local coordinates with the maximum range on each 
design variable of -1 to + 1. The program then 

4. asks the user to enter an integer which specifies whether design point coordinates are 
to be also generated in global coordinates. If they are to be calculated in global 
coordinates, the program then 

5. prompts the user to enter the range of design variables in global coordinates. 

Results with commentary are written to file "design.res". Design points without commentary 
are written to file "design.run". 

8J Augmented Minimum Point Design for 3rd and 4th Order Approximation 

A 3 k factorial design is used as the building block of this minimum point design. The 3 k 
factorial design provides information for calculating the coefficients associated with the 
terms 

1, x p xpc., xf, x]x p xfxj, j*i (104) 

Additional points are then added at -1 and 1 (in local coordinates) along the x, axis. These 
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points together with the 3 k factorial design point give the star pattern which can be seen in 
Figure 31. With this arrangement of points, there are 5 design points along the x* axis which 
provides information for calculating the coefficient associated with the terms 

< ( 105 ) 


Additional design points are then placed in each Xj, Xj plane which provides information for 
calculating the coefficient associated with the terms 


3 

X:X. 


( 106 ) 


These points are also shown in Figure 31. 

8.3.1 Specifics of program DESIGN4 

A listing of the FORTRAN program DESIGN4 is found in Appendix 2 and a copy of that 
program is found in file "design4.f on the floppy disk accompanying this report. The 
program should be compiled with a F77 compiler with the compiled program called 
"design4". To run the program just enter "design4 H from the keyboard. The program 
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prompts the user for needed information. Prompts and response are similar to those for the 
program DESIGNS. 


8.4 Conclusion 

A minimum point design is a design that has just enough design points to allow the 
determination of the coefficients of an approximating polynomial. An augmented minimum 
point design is a minimum point design which contains additional design points. Augmented 
minimum point designs are competitive with, if not better than, central composite designs 
for developing a 2nd order response surface. Minimum point designs should be augmented 
with enough points that the approximation is 20-50% over-determined. A program 
DESIGNS was developed for generating augmented minimum point designs for developing 
a 2nd order response surface. 

When there are more than 6 design variables, 3rd or higher order approximations require 
so many design points that it is computationally better to perform a sequence of 2nd order 
approximations in an optimization process than one higher order approximation. When 
there are 6 or fewer design variables, a 2nd order approximation may often be satisfactory. 
However, for those cases where it is desirable to use a higher order approximation, program 
DESIGN4 was developed. That program generates designs which can be used to develop 
3rd or 4th order approximations. 
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9. Solution Algorithm 


9.1 Introduction 

In this investigation, the program NEWPSI was used to perform the studies involving 
polynomial approximations. That program can investigate under-determined, exactly- 
determined, or over-determined approximations of various orders. The version submitted 
with this report can handle up to 15 design variables as programmed. The order of 
polynomial it can handle is as follows: 

1. one design variable, up to a 20th order polynomial 

2. two design variables, up to a 5th order polynomial 

3. for 2-15 design variables, a second order polynomials. 

One can use up to 250 designs to train the approximation. In calculating v G , it can handle 
up to 2000 grid points. 

The program solves for the undetermined parameters associated with the approximation. 
It then evaluates the approximate function at the design points and calculates the error 
parameter, v. It then reads in the design points and function value on the test grid. The 
approximate function is evaluated at the grid points and the error parameter, v G , is then 
evaluated. 

9.2 Program Specifics 

A listing of the FORTRAN program NEWPSI is found in Appendix 3 and a copy of that 
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program is found in file "newpsi.f' on the floppy disk accompanying this report. The 
program should be compiled with a F77 compiler and the compiled program called "newpsi". 
To run the program just enter "newpsi" from the keyboard. Data is read from the file 
"newpsi.dat". Data can be in free format. The program asks for the following data: 

1. a value of the print code, ip; (If ip = 4, great quantities of output are generated for 
program checkout. Normally the program is run with ip=0 for normal output). 

2. the number of design variable, nd; 

3. the order of the polynomial being considered, np; 

4. the number of design points in the design, m; 

5. the design and function value at the design points, x(ij), y(i); 

6. the number of design points on the grid, ng; and 

7. the design and function value at the grid points, xx(i,j), yy(j). 

Output is written to the screen and to file "newpsi.res". 
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10. Conclusion 


For a given order of approximation of a function, f, the quality of the approximation is 
affected by 

a. the number of levels of the design variables, 

b. the location of the design points, and 

c. the degree which the approximation is over-determined. 


For an nth order approximation, 

1. there must be n + 1 levels of the design variables; 

2. the design points must be located so that information is available for calculating all 
of the nth derivatives of f; 

3. the approximation should be, at least, 20-50% over-determined. 


For example, for a 2nd order approximation in 3 design variables, there must be at least 3 
levels of the design variables, design points must be located so that information is available 
for calculating 


M. JL, 

dx’ dxifo) 


<■1,3; 7=1,3 


(107) 


A complete 2nd order polynomial approximation contains 10 undetermined coefficients. 
Thus, at least 10 design points are required to provide information for calculating these 
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coefficients. To have the approximation 30% over-determined, one would want to use 13 
design points. 

For second order approximations, when there are fewer than 6 design variables, central 
composite designs meet requirements 1-3. However, for 6 or more design variables, these 
designs contain too many design points. A minimum point design is one which contains just 
enough design points, meeting the derivative requirements of item 1 and 2 above, to exactly- 
determine the approximation. An augmented minimum point design is a minimum point 
design supplemented with additional design points. The program DESIGNS was developed 
to yield augmented minimum point designs for 2nd order approximations. The quality of 
approximations developed using designs from program DESIGNS was comparable to, if not 
better than, other standard designs such as the central composite designs. 

For more than 6 design variables, 3rd and 4th order approximations require so many design 
points to determine the coefficients in those approximations that it is more computationally 
efficient to develop a number of 2nd order approximations than one approximation of 3rd 
or higher order. For 6 or fewer design points, 2nd order approximations may be quite 
adequate. However, for those cases where one wishes to use a 3rd or 4th order 
approximation, the program DESIGN4 was develop. That program generates an augmented 
minimum point design for developing a 4th order approximation. 


Previous studies have shown that the quality of approximations using neural networks is 
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comparable to those using polynomial approximations when the number of undetermined 
parameters associated with the approximations is the same. Thus, neural networks trained 
with designs from DESIGNS or DESIGN4 should offer approximations of comparable 
quality to those obtained using polynomial approximations with the same number of 
associated undetermined parameters. 
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Figure 1. Artificial neural net 
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Figure 2. Complete 
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Figure 3. Deficient design 
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Figure 6. Design from program DESIGNS 
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Figure 8. Fox 




One Dimensional Example 
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Figure 9. One dimensional example 
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Figure 10. D, A, E, and G optimality, first order approximation 
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Parameter vG 

Second Order Polynomial Approximation 



— Uniform — 
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-A— E Optimality — * 

G Optimality 


Figure 11. D, A, E, and G optimality, second order approximation 
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Figure 12. D, A, E, and G optimality, third order approximation 
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Parameter vG 


Parameter vG 

Fourth Order Polynomial Approximation 
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Figure 13. D, A, E, and G optimality , fourth order approximation 
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Parameter vG 


Parameter vG 

First Order Polynomial Approximation 
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Parameter vG 

Second Order Polynomial Approximation 



Uniform — t— S-Optimality Q Optimality 


Figure 15. S and Q optimality, second order approximation 
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Parameter vG 

Third Order Polynomial Approximation 
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Figure 16 . 


S and Q optimality, third order approximation 
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Parameter vG 


Parameter vG 

Fourth Order Polynomial Approximation 



Figure 17. S and Q optimality, fourth order approximation 
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Y and its Approximation 

Q optimality, 1 1 points out of 13 
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Figure 18. Q optimality, 11 out of 13 points selected 
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Figure 19. Q optimality, 7 out of 13 points selected 





Y and its Approximation 

Q optimality, 5 out of 13 points 
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igure 20. Q optimality, 5 out of 13 points selected 
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First Order Polynomial Approximation 



Random points, first order approximation 
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Figure 21. 



Parameter vG 


Parameter vG 

Second Order Polynomial Approximation 



Figure 22. Random points, 


second order approximation 
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Parameter vG 


Parameter vG 

Third Order Polynomial Approximation 



Figure 23. Random points, third order approximation 
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Parameter vG 


Parameter vG 

Fourth Order Polynomial Approximation 



Figure 24. Random points, fourth order approximation 
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MEASUREMENTS OF QUALITY OF FIT 
BEFORE AND AFTER t-test 
PERFORMED ON 

y = 10 x[ - 2 0X2 x\ + IOX 2 X 1 +xf — 2xi + 5 


"Fox's Banana Function" 

SECOND ORDER APPROXIMATION 

Y= b 0 + Jbi xi + b 2 x z + b 2 xl + i><x ix 2 + £ 5*2 

Before t-test After t-test 


v : 
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v ; 
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Figure 25. Significance testing, Example 1, 2nd order 
approximation 
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MEASUREMENTS OF QUALITY OF FIT 
BEFORE AND AFTER t-test 
PERFORMED ON 

y= 10x* - 20x2x1 + 10 x\x\+ x\ -2xi + 5 
"Fox's Banana Function" 


THIRD ORDER APPROXIMATION 

Y — bo + bi Xi + biX2 + jb jX[ + £>4X1X2 + £5X7 + jbgx^ + b~i x 1 X2 + bn x 1x5 + bg X2 
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Figure 26. 


Significance testing, Example 1 
approximation ' 


3rd order 
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measurements of quality of fit 

BEFORE AND AFTER t-test 
PERFORMED ON 

Y = (4 + Xi ) 3 + sin[^*(x! + l)] + 2 + x^ + sinjj-) + 7x 2 x, 


SECOND ORDER APPROXIMATION 

Y= b 0 + Jbi X! + b 2 x 2 + b 3 xf + D4X1X2 + b 5 xl 


Before t-test 


After t-test 
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Figure 27 . 


Significance testing, Example 2, 

approximation 


2nd order 
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MEASUREMENTS OF QUALITY OF FIT 
BEFORE AND AFTER t-test 
PERFORMED ON 


Y = (4 + Xi) 3 4- sin 


-y * (xi + 1 


>] 


4- 2 + xi 4* sin 


(f) 


+ 7x 2 Xi 


THIRD ORDER APPROXIMATION 

Y- b 0 + bi xi + bzXz + b 3 Xi + i> 4 xix 2 + b$x\ + b&xf + bixlx 2 + baXixj + b 9 x| 


Before t-test 


After t-test 


v : 0.7 
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Figure 28. Significance testing, Example 2, 3rd order 
approximation 
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Figure 29. The 35 bar truss 
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Figure 30. Additional points to complete a second order design 



added to find coefficients of terms of Eq . (1 06) 
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Figure 31. Additional points to complete a fourth order des g 
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Appendix 1 

Program DESIGNS 

PROGRAM DESIGNS 
C 

C PROGRAM TO GENERATE DESIGNS FOR 2ND ORDER POLYNOMIAL 

C PROGRAM DIMENSIONED FOR UP TO 20 VARIABLES 

C RESULTS TO SCREEN AND TO FILE designs. res 

C DESIGN IN GLOBAL COORDINATES TO FILE designs. run 

C 

C DEFINITIONS 

C N = NUMBER OF DESIGN VARIABLES 

CM = NUMBER OF RANDOM DESIGNS POINTS 

C 

DIMENSION X(2000, 20) 

DIMENSION XBB (20) ,XBE(20) ,A(20) ,B(20) 

1 FORMAT (I5,6F10.6) 

2 FORMAT ( ' PROGRAM GENERATES DESIGNS FOR FITTING 2ND ORDER', 

X' POLYNOMIAL') 

3 FORMAT (' ENTER NUMBER OF DESIGN VARIABLES') 

4 FORMAT ( ' NUMBER OF DESIGN VARIABLES = N =' , 13) 

11 FORMAT (6F10. 6) 

OPEN (UNIT=7 , FILE=' designs . res ' ) 

OPEN (UNIT=8 , FILE= ' designs . run ' ) 

WRITE (6, 2) 

WRITE (6, 3) 

READ ( 5 , * ) N 
WRITE (6, 4) N 
C SET UP TERMS 

NP1=N+1 
NM1=N-1 

M=(N*N+3*N+2) / 2 
MP1=M+1 

ZERO DESIGN MATRIX 
DO100I=l,M 
DOIOO J=1 , N 

100 X(I,J)=0. 

11=0 

c 

C * ’ * 

C GENERATE THE FIRST N+l POINTS FOR FITTING A LINEAR FUNCTION 

C THE FIRST POINT IS WHEN ALL X'S ZERO, ALREADY DONE 

C GENERATE NEXT N POINTS 

D0101I=1 , N 
11=1+1 

101 X(II,I)=1. 

C 

c 

C GENERATE NEXT N POINTS 

C THE 2N+1 POINTS THUS GENERATED WILL ALLOW ADDING SQUARED TERMS 
C 

DO102I=l , N 
II=I+N+1 

102 X(II,I)=-1. 

C 

C 

c 

C GENERATE NEXT N(N-l)/2 POINTS 

C THE (N*N+3 *N+2 ) / 2 POINTS THUS GENERATED WILL ALOW ADDING CROSS 

C PRODUCT TERMS. WE WILL THEN HAVE COMPLETE 2ND ORDER POLYNOMIAL 

C APPROXIMATION 

C 


123 



on o n o o n 


ILAST=2*N+1 
IDO=N- 1 


J=1 

JJ=2 

103 CONTINUE 
DO104I=l, IDO 
II=I+ILAST 
X(II,J)=1. 
X(II,JJ)=1- 
JJ=JJ+1 

104 CONTINUE 
ILAST=ILAST+IDO 
IDO=IDO- 1 
J=J+1 


C 

C 

C 


C 

C 

C 


JJ=J+1 

IF ( J . LE . NM1) GOTO103 


IF WE GOT HERE WE HAVE DEVELOPED THE MINIMUM POINT DESIGN 

SS3{?-:}: s »£ S £SS 5 ™ £5 £ SS3Z 

WRITE ( 6 | * ) ' DESIGN POINTS WRITTEN TO FILE designs. res' 

nFVELOP DESIGN POINTS TO AUGMENT THE MINIMUM POINT DESIGN 
THE NUMBER OF RANDOM DESIGN POINTS TO BE DEVELOPED 
WRITE (6 , *) ' ENTER THE NUMBER OF RANDOM GENERATED DESIGN PT , 

X' DESIRED=M' 

READ ( 5 * ) M 

WRITE (6 *)' NUMBER OF RANDOM DESIGN POINTS M=' ,M 

WRITE d' *) ' NUMBER OF RANDOM DESIGN POINTS M— ' ,M , 

IFLAG IS ANY POSITIVE INTEGER USED TO START RANDOM', 

X' PROCESS' 

WRITE (6,*)' ENTER IFLAG' 

READ ( 5 , * ) IFLAG 
WRITE (6,*)' IFLAG= ' , IFLAG 
WRITE (7,*)' IFLAG= ' , IFLAG 
DO850I=l,M 
11 = 11+1 


D0851J=1,N 
I FLAG= I FLAG + 1 
XDUM=RAND ( I FLAG ) 
X(II, J)=2.*XDUM-1. 
851 CONTINUE 
850 CONTINUE 


IF WE GOT HERE WE HAVE FINISHED 
WRITE(6,*)' RANDOM DESIGN POINTS 


GENERATING THE RANDOM DESIGN PTS 
WRITTEN TO FILE designs. res' 


PRINT OUT THE MINIMUM POINT MATRIX IN LOCAL COORDINATES 


700 


WRITE ( 7 , * ) ' DESIGN MATRIX IN LOCAL COORDINATES' 
ITOTAL=II 


D07 001 = 1 , ITOTAL 

WRITE (7,1)1, (X(I,J) , J=1 , N) 

SEE IF WE ARE TO GENERATE DESIGNS IN GLOBAL COORDINATES 


WRITE ( 6 , * ) ' ITEST=1 IF DESIGN POINTS ARE TO BE IN GLOBAL', 
X' COORDINATES' 

WRITE(6,*)' OTHERWISE, ITEST=0' 
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WRITE (6,*)' ENTER ITEST' 

READ (5, *) ITEST 
IF (ITEST . NE. 1 ) GOTO860 

C IF WE GOT HERE WE ARE TO GENERATE DESIGNS IN GLOBAL COORDINATES 

WRITE ( 6 , * ) ' ENTER LOWER AND UPPER RANGE ON EACH DESIGN VARIABLE' 
WRITE (6,*)' i.e. ENTER XBB(I) TO XBE(I)' 

D0861I=1 , N 

READ ( 5 , * ) XBB (I) , XBE ( I ) 

WRITE (6,*) ' I, XBB (I) ,XBE(I)=' , I, XBB (I) , XBE(I) 

WRITE ( 7 , * ) ' I , XBB (I) , XBE ( I ) = ' , I , XBB (I) , XBE ( I ) 

861 CONTINUE 
G0T0862 

860 CONTINUE 

IF WE GOT HERE LOWER BOUND VARIABLE IN GLOBAL COORDINATES IS -1 
IF WE GOT HERE UPPER BOUND VARIABLE IN GLOBAL COORDINATES IS 1 
D0863I=1 , N 
XBB ( I ) =-l . 

XBE ( I ) =1 . 

863 CONTINUE 

862 CONTINUE 

WRITE ( 7 , * ) ' I , XBB (I) , XBE (I),A(I),B(I)' 

D013 011=1 , N 

A ( I ) = (XBE ( I ) -XBB ( I ) ) / 2 . 

B ( I ) = (XBE ( I ) +XBB ( I ) ) 12 . 

WRITE ( 7 , * ) I , XBB ( I ) ,XBE(I) ,A(I) ,B(I) 

1301 CONTINUE 

D01202I=1 , ITOTAL 
DO1202 J=1 , N 

1202 X(I, J)=A(J) *X(I, J)+B(J) 

WRITE ( 6 , * ) ' DESIGN IN GLOBAL COORDINATES WRITEN TO designs. res' 
WRITE ( 6 , * ) ' DESIGN IN GLOBAL COORDINATES WRITEN TO designs. run' 
WRITE (7,*)' DESIGN IN GLOBAL COORDINATES' 

WRITE ( 8 , * ) ITOTAL 

D0970I=1 , ITOTAL 

WRITE (7,1)1, (X(I, J) , J=1 , N) 

WRITE(8 , 11) (X(I,J) , J=1 , N) 

970 CONTINUE 
STOP 
END 
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PROGRAM DESIGN4 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 

C 


PROGRAM TO GENERATE DESIGNS FOR 4TH ORDER POLYNOMIAL 
PROGRAM DIMENSIONED FOR UP TO 6 VARIABLES 
RESULTS TO SCREEN AND TO FILE design4*:res 
DESIGN IN GLOBAL COORDINATES TO FILE design4.run 


DEFINITIONS 

N = NUMBER OF DESIGN VARIABLES 

M = NUMBER OF RANDOM DESIGNS POINTS 


4 

11 


DIMENSION X (2000 , 6 ) 

DIMENSION XBB(IO) ,XBE(10) ,A(10) ,B(10) 

1 FORMAT (15 , 6F10 . 6) 

2 FORMAT (' PROGRAM GENERATES DESIGNS FOR FITTING 4TH ORDER , 
X' POLYNOMIAL') 

3 FORMAT (' ENTER NUMBER OF DESIGN VARIABLES') 

FORMAT ( ' NUMBER OF DESIGN VARIABLES = N =' , 13) 

FORMAT (6F10. 6) 

OPEN (UNIT=7 , FILE= ' design4 . res ' ) 

OPEN (UNIT=8 , FILE= ' design4 . run ' ) 

WRITE (6, 2) 

WRITE (6,3) 

READ ( 5 , * ) N 
WRITE ( 6 , 4 ) N 
IF (N . EQ . 6 ) GOTO601 
IF (N . EQ . 5 ) GOTO501 
IF (N . EQ . 4 ) G0T04 0 1 
IF (N . EQ . 3 ) G0T03 01 
IF (N . EQ . 2 ) G0T02 0 1 
IFfN.EQ. 1) GOTOlOl 

WRITE (6,*)' PROGRAM CAN NOT DO MORE THAN 6 DESIGN VARIABLES 
WRITE (7 , * ) ' PROGRAM CAN NOT DO MORE THAN 6 DESIGN VARIABLES 
STOP 


t 

t 


DEVELOP 3 FACTORIAL DESIGN TO GET 4 DESIGN VARIABLE PRODUCT TERMS 
101 CONTINUE 
11=0 

D0100I1=1, 101, 50 
11 = 11+1 

X (II , 1) =FLOAT (11-51) /100. 

100 CONTINUE 
GOTO701 
201 CONTINUE 
11=0 

D0200I1=1, 101, 50 
D0200I2=1, 101, 50 
11 = 11+1 

X ( II , 1 ) =FLOAT (11-51) / 100. 

X ( II , 2 ) =FLOAT (12-51) /100. 

200 CONTINUE 
GOTO701 
301 CONTINUE 
11=0 

D0300I1=1, 101, 50 
D0300I2=1 , 101,50 
D0300I3=1 , 101,50 
11 = 11+1 

X(II, 1) =FLOAT (11-51) / 100. 

X(II, 2)=FLOAT(I2-51) /100. 
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non 


X(II, 3)=FLOAT(I3-51) /100. 

300 CONTINUE 
GOTO701 
401 CONTINUE 
11=0 

D0400I1=1, 101, 50 
D0400I2=1, 101, 50 
D0400I3=1, 101, 50 
D0400I4=1, 101, 50 
11=11+1 

X(II, l)=FLOAT( 11-51) /100. 

X (II , 2 ) =FLOAT (12-51) /100. 

X ( I I , 3 ) =FLO AT (13-51) /100. 

X(II , 4 ) =FLOAT (14-51) /100. 

400 CONTINUE 
GOTO701 
501 CONTINUE 
11=0 

D0500I1=1, 101, 50 
D0500I2=1, 101, 50 
D0500I3=1, 101, 50 
D0500I4=1 , 101,50 
D0500I5=1, 101, 50 
11=11+1 

X(II, l)=FLOAT (11-51) /100. 

X(II,2)=FLOAT (12-51) /100. 

X ( II , 3 ) =FLOAT ( 13-51 ) /100. 

X ( I I , 4 ) =FLOAT (14-51) /100. 

X(II,5)=FL0AT( 15-51) /100. 

500 CONTINUE 
GOT0701 

601 CONTINUE 
11=0 

D0600I1=1 , 101, 50 
DO600I2=l, 101, 50 
D0600I3 = 1, 101, 50 
D0600I4=1 , 101, 50 
D0600I5=1 , 101,50 
D0600I6=1, 101, 50 
11=11+1 

X(II,1)=FL0AT (11-51) /100. 

X ( II , 2 ) =FLOAT (12-51) /100. 

X (II , 3 ) =FLOAT (13-51) /100. 

X (II , 4 ) =FLOAT ( 14 -51 ) /100. 

X (II , 5 ) =FLOAT (15-51) /100. 

X ( II , 6 ) =FLOAT (16-51) /100. 

600 CONTINUE 
G0T07 0 1 

701 CONTINUE 

ENTER REST OF POINTS IN THE STAR FORMATION 

D0702I=1 , N 
11=11+1 
DO703 J=1 , N 
703 X(II,J)=0. 

X(II,I)=1. 

702 CONTINUE 
DO704I=l , N 
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o o o n 


11 = 11+1 
DO705 J=1 , N 
705 X(II,J)=0. 

x(ii,i)=-i. 

704 CONTINUE 

enter terms to calculate coefficient associated with the term 
X ( I ) **3*X ( J) 

NM1=N-1 

IDO=N-l 

J=1 

JJ=2 

803 CONTINUE 
D0804I=1 , IDO 
11 = 11+1 
X(II, J)-l. 

X(II, JJ)=-5 
11 = 11+1 
X(II, J)=-5 
X(II, JJ)=1. 

JJ=JJ+1 

804 CONTINUE 
IDO=IDO-l 
J=J+1 
JJ=J+1 

IF ( J . LE . NM1) GOTO803 


C 

C 

c 


c 

c 

c 


c 

c 

c 

c 


IF WE GOT HERE WE HAVE DEVELOPED THE MINIMUM POINT DESIGN 

S™'*.’,: w w e e £3 S3SS :;S:: 333 3 3! S 3 K 

WRITE DESIGN POINTS WRITTEN TO FILE design4.res' 

DFVELOP DESIGN POINTS TO AUGMENT THE MINIMUM POINT DESIGN 
SSn ?N THE NUMBER OF RANDOM DESIGN POINTS TO BE DEVELOPED 
WRITE ( 6 , * ) ' ENTER THE NUMBER OF RANDOM GENERATED DESIGN PTS', 

X' DESIRED=M' 

READ ( 5 *)M 

WRITE ( 6 *) ' NUMBER OF RANDOM DESIGN POINTS M =/ ,M 

WRTTF n' ' NUMBER OF RANDOM DESIGN POINTS M— ' ,M V mAM / 

Sme(6:*)' “ is aot positive integer used to start random' , 

X' PROCESS' 

WRITE (6,*)' ENTER IFLAG' 

READ (5, *) IFLAG 
WRITE (6,*)' IFLAG=' , IFLAG 
WRITE (7,*)' IFLAG= ' , IFLAG 
DO850I=l,M 
11=11+1 
D0851J=1 , N 
IFLAG=IFLAG+1 
XDUM=RAND ( I FLAG ) 

X(II , J) =2 . *XDUM-1 . 

851 CONTINUE 
850 CONTINUE 

IF WE GOT HERE WE HAVE FINISHED GENERATING THE RANDOM DESIGN PTS 
WRITE ( 6^* ) 'RANDOM DESIGN POINTS WRITTEN TO FILE design4 . res ' 

PRINT OUT THE MINIMUM POINT MATRIX IN LOCAL COORDINATES 
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n o o 


c 


WRITE (7 , *) ' DESIGN MATRIX IN LOCAL COORDINATES' 
ITOTAL=II 


DO700I=l, ITOTAL 
WRITE (7,1)1, (X(I, J) , J=1 , N) 
700 CONTINUE 


C 

C 


SEE IF WE ARE TO GENERATE DESIGNS IN GLOBAL COORDINATES 

WRITE (6 , *) ' ITEST=1 IF DESIGN POINTS ARE TO BE IN GLOBAL' 

X' COORDINATES' ' 

WRITE ( 6 , * ) ' OTHERWISE, ITEST=0' 

WRITE (6,*)' ENTER ITEST' 

READ ( 5 , * ) ITEST 

IF ( ITEST. NE. 1) GOTO860 

WRTTWfi°*\ ARE T0 GENERATE DESIGNS IN GLOBAL COORDINATES 

, ? NTER L0WER ^ UPPER RANGE ON EACH DESIGN VARIABLE' 
WRITE (6,*)' i.e. ENTER XBB(I) TO XBE (I) ' 

D0861I=1 , N 

READ ( 5 , * ) XBB (I) , XBE ( I ) 

WRITE (6, *) ' I,XBB(I) , XBE (I ) -',1 , XBB (I) ,XBE(I) 

861 otJmnUe* 1 ” I ' XBB(I) ' XBE < I ) = '' I ' XBB < i ).™E(I) 

GOT0862 
860 CONTINUE 


LOWER 

UPPER 


BOUND 

BOUND 


VARIABLE 

VARIABLE 


IN 

IN 


I , XBB (I) , XBE ( I ) , A (I ) , B ( I ) 


IF WE GOT HERE 
IF WE GOT HERE 
D0863 1=1 , N 
XBB ( I ) =-l . 

XBE (I) =1 . 

863 CONTINUE 
862 CONTINUE 

WRITE (7 , *) ' 

DO1301I=l , N 
A (I) = (XBE (I) -XBB (I) )/2. 
B(I)=(XBE(I)+XBB(I) ) / 2. 

WRTTE ( 7 , * ) 1 , XBB ( I ) , XBE ( I ) ,A(I) ,B(I) 

1301 CONTINUE 1 ' 

D01202I=1, ITOTAL 
DO1202 J=1 , N 

X ( I , J) =A ( J) *X ( I , J) +B ( J) 

WRITE (6,*)' DESIGN IN GLOBAL COORDINATES 

WRITE ( 6 , * ) ' DESIGN IN GLOBAL COORDINATES 

WRITE (7 , * ) ' DESIGN IN GLOBAL COORDINATES 

WRITE (8 ,*) ITOTAL 

D0970I=1, ITOTAL 

WRITE (7,1)1, (X (I , J) , J=1 , N) 

WRITE (8, 11) (X (I , J) , J=1 , N) 

970 CONTINUE 
STOP 
END 


GLOBAL COORDINATES IS -1 
GLOBAL COORDINATES IS 1 


1202 


WRITEN TO desigr>4 . res ' 
WRITEN TO design4.run' 
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Appendix 3 
program NEWPSI 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM newpsi 


Z STSS^SSS: ex^arof oei? r ae^inea UhiC ^ 

it^can hanaie up to 15 aesign a = 

r sr o„f asssr^.^^si*" 

one can use up to’lSO aesigns to train the approximation. 

It can handle up to 2000 grid points 

IMPLICIT REAL* 8 (A-H,0-Z) 

dimension x(250,15) ,y(250) ,a(250,136) 

dimension yhat(250) 

dimension xx(2000,15) ,yy(2000) , abig (2000 , 136) 
dimension yyhat(2000) 

. FORMAT (9F8 . 4) 

; FORMAT ( Flo? 6^ 1H , ,P10.6,1H,/F10.6 / 1H,, F10. 6, 1H, , F10. 6, 1H, ,F10 

X1H, , F10 . 6) . ... 

OPEN (UNIT=5 , FILE= ' newpsi - dat ) 

OPEN (UNIT=7 , FILE= ' newpsi . res ' ) 

OPEN (UNIT=8,FILE=' newpsi. plot') 

*************************************************************** 


read in data 


read in the print code 
read (5, *) ip 

enter number of design variables, nd 
read ( 5 , *) nd 

enter THE DEGREE OF POLYNOMIAL TO BE CONSIDERED, np 
READ ( 5 , * ) np 

ENTER NUMBER OF DESIGNS FOR PROBLEM, M 
READ ( 5 , * ) M 


write ( 6 , *) ' 
write(6, *) ' 


print code ip-'»ip 

number of design variables, nd- ,na 
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write(6,*)' degree of polynomial being considered=np=' , np 
write (6 ,*)' number of designs m=' ,m 
write (7,*)' print code ip=',ip 

write(7,*)' number of design variables, nd=',nd 
write(7,*)' degree of polynomial being considered=np=' , np 
write(7,*)' number of designs m=',m 
C 

c read in designs and set up matrix a 

C 

write (7 , *) ' x(i, j) ,y(i) ' 

D0101I=1,M 

read (5 , *) (x(i, j) , j=l,nd) ,y(i) 
write (7,*) (x(i, j) ,j=l,nd) ,y(i) 

101 continue 
c 

c set up the coefficient matrix, a, in the matrix equation 

c y=a x 

c 

call geta (ip,m, nd, np, n, x, a) 

C 

C SEE WHETHER SYSTEM IS UNDER, EXACTLY , OR OVER DETERMINED 

C 

IF (M. GE . N) GOTO400 

C IF WE GOT HERE WE ARE UNDER-DETERMINED 

WRITE (6,*)' SYSTEM IS UNDER-DETERMINED' 

WRITE ( 7 , * ) ' SYSTEM IS UNDER-DETERMINED' 

CALL PSI (ip,M,N,A,Y,B) 

GOTO402 

400 CONTINUE 
IF(M.GT.N) GOTO401 

C IF WE GOT HERE WE ARE EXACTLY DETERMINED 

WRITE ( 6 , * ) ' SYSTEM IS EXACTLY DETERMINED' 

WRITE (7 , *) ' SYSTEM IS EXACTLY DETERMINED' 

CALL EXACT (ip,M,A,Y,B) 

GOTO402 

401 CONTINUE 

C IF WE GOT HERE WE ARE OVER-DETERMINED 

WRITE (6,*)' SYSTEM IS OVER-DETERMINED' 

WRITE ( 7 , * ) ' SYSTEM IS OVER-DETERMINED' 

CALL OVER (ip,M,N,A,Y,B) 

402 CONTINUE 
C 

C 

C EVALUATE APPROXIMATION AT DESIGNS 

C 

WRITE (6,*)' MATRIX OF COEFFICIENTS, B(I)' 

WRITE (7,*)' MATRIX OF COEFFICIENTS, B(I)' 

WRITE(6,*) (B(I) ,I=1,N) 

WRITE ( 7 , * ) (B ( I) ,1=1, N) 

WRITE (7,*)' MATRICES Y ( I ) AND YHAT(I) ' 

C 

c recalculate matrix a 

call geta (ip,m, nd, np, n, x, a) 
c 

c calculate approximation at designs and print results 

c 

write ( 7 , * ) ' y ( i ) , yhat ( i ) ' 

D0102 1=1 , M 
YHAT ( I ) =0 . 

DO103 J=1 , N 
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o o 


103 

102 


601 


603 


c 

c 

c 

602 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


yhat ( i ) =yhat (i)+a(i,j)*b(j) 

CONTINUE 

WRITE (7 , *) Y (I) , YHAT ( I ) 

CONTINUE 

evaluate function at grid 
read(5,*)ng 

write(6,*)' number of designs on grid = ngn',ng 
write(7,*)' number of designs on grid = ngn',ng 
write (7 , *) ' xx ( i , j ) , yy ( i ) ' 

DO601I=l , ng 

read (5,*) (xx ( i , j ) , j=l , nd) , yy ( i) 
write (7 , *) (xx(i, j) , j=l,nd) ,yy(i) 
continue 

call getabg ( ip , ng , nd , np , n , xx , abig ) 
write (7 , *) ' yy (i) /YY^at (i) at grid' 

D0602I=1 , ng 
YYHAT ( I ) =0 . 

DO603 J=1 , N 

yyhat ( i ) ^yyhat ( i ) +abig ( i , j ) *b ( j ) 

CONTINUE 

WRITE ( 7 , * ) YY ( I ) , YYHAT (I) 

write the plot file 

write (8 , *) (xx (i, j ) , j=l/ nd) , yyhat (i) 

CONTINUE 

calculate statistical terms 

call stat it (m , y , yhat , ng , yy , yyhat) 

STOP 

END 

subroutine geta ( ip , m , nd , np , n , x , a ) 

************************************************************* 

************************************************************* 

This subroutine generates the matrix a where the matrix 
equation is y= a b. Here y are the training functions, 
b are undetermined coefficients. The algorithm is programmed 
to handle 

1. any level of approximation for one design variable 

2. up to 5th order polynomial in two design variables 

3. quadratic approximation in more than two design variabaales 

************************************************************ 

************************************************************ 


IMPLICIT REAL* 8 (A-H,0-Z) 
dimension x(250, 15) ,a(250, 136) 

do for each design 

do300i=l,m 

******************************************************** 
if nd is not equal to 1 go to 400 
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if (nd.ne. I)goto400 
c 

q ★*********************************ifc-****TkTlr******* , * , *-*-* , *'* , **r* 

C ************************************ *-* ****************** 

C 

c here we have nd=l, i.e. one design variable 

c we will develp a's for all np's 

c 

a ( i , 1) =1 • 
j=l 

do201k=l, np 
j=j+l 

a(i, j)=x(i,l)**k 
201 continue 
n=np+l 
goto301 
c 
c 

C ★★★★★★★★★★★★★★★★’A******************* 

400 continue 
c 

c if nd is not equal to 2 go to 500 

if (nd . ne • 2 ) goto500 
c 



c 

c if we got here we have 2 design variables 

c 

xl=x(i, 1) 
x2=x(i, 2) 
c 

c ********************* 

c 

c add the constant and linear terms 

c 

a (i , 1) =1 . 

a (i, 2 ) =xl 
a (i, 3) =x2 
n=3 

if (np. It . 2) got o301 
c 

c ******************** 

c 

c add the 2nd order terms 

c 

a (i , 4 ) =xl**2 
a (i , 5 ) =xl*x2 
a( i , 6) =x2**2 
n=6 

if (np. It. 3) goto301 
c 

c ******************* 

c 

c add the cubic terms 

c 

a ( i , 7 ) =xl**3 
a ( i , 8 ) =xl**2*x2 
a (i , 9 ) =xl*x2**2 
a(i, 10)=x2**3 
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n=10 

if (np. It. 4) goto301 
c 

c ******************* 

c 

c add the 4th order terras 

c 

a(i, 11) =xl**4 
a(i,12)=xl**3*x2 
a(i, 13 ) =xl**2*x2**2 
a(i, 14)=xl*x2**3 
a (i , 15) =x2**4 
n=15 

if (np. It. 5)goto301 
c 

c ******************* 

c 

c add the 5th order terras 

c 

a (i, 16) =xl**5 
a(i, 17)=xl**4*x2 
a(i, 18) =xl**3*x2**2 
a ( i , 19 ) =xl**2*x2**3 
a(i,20)=xl*x2**4 
a ( i , 21) =x2**5 
n=21 

if (np. It, 6) goto301 
c 

c ****************** 

c 

c algorithm not programed for polynomials of order larger than 5 

c 

write(6,*)' for two design variables, algorithm not programed for' 
write(6,*)' polynomials of order larger than 5' 

write(7,*)' for two design variables, algorithm not programed for' 
write(7,*)' polynomials of order larger than 5' 
stop 
c 

c ****************************************************** 

c ****************************************************** 

c 

500 continue 
c 

c if we got here number of design variables >2 

c 

c ******************** 

c 

c enter constant and linear terms 

c 

a (i, 1)=1* 
j=l 

do501k=l , nd 
j=j+l 

a(i, j)=x(i,k) 

501 continue 

n=j 

if (np. It.2)goto301 
c 

Q ******************** 

c 
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c 

c 


enter the quadratic terms 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


do502k=l , nd 
do 5 02L=k f nd 
j=j+l 

a (i f j)=x(i,k) *x(i,L) 
502 continue 

n=j 

if ( np . 1 1 . 3 ) goto3 0 1 
******************** 


algorithm not programmed for more than quadratic approximation 
when number of design variables >2 

wr f ^ e (®»*) ' algorithm not programmed for more than quadratic' 

' approximation when number of design variables >2' 
write (7 , *) ' algorithm not programmed for more than quadratic' 
*) ' approximation when number of desiqn variables >2' 

stop 




**************** 

**************** 


********* 

******** 


print out some results 


301 continue 

if (ip.lt. 4)goto302 

write (6, *) ' a (i, j ) ', (a ( i , j ) , j=l , n) 

write (6 , *) ' ' 

write (7 , *) ' a ( i , j ) (a (i , j ) , j=i , n) 
write (7 , *) ' ' 

302 continue 


c 

c 

c 






c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


300 continue 

wr fte(6,*)' number of undetermined coef=n=' / n 
write (7,*)' number of undetermined coef=n=',n 

return 

end 

subroutine getabg ( ip , m , nd , np , n , x , a ) 






This subroutine generates the matrix a where the matrix 
equation is y= a b. Here y are the training functions, 

to a handle etermined coeff icients * The algorithm is programmed 


1. any level of approximation for one design variable 

2. up to 5th order polynomial in two design variables 

3. quadratic approximation in more than two design variabaales 


*************** 

*************** 




IMPLICIT REAL*8 (A-H,0-Z) 
DIMENSION A (2000 ,136) 
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DIMENSION X ( 2 000 , 15 ) 


c 

c do for each design 

c 

do300i=l,m 

c 

c ******************************************************** 

c 

c if nd is not equal to 

if (nd.ne, l)goto400 
c 

c ********************* 

c ********************* 

c 

c here we have nd=l, i.e. one design variable 

c we will develp a's for all np's 

c 

a ( i , 1) =1 • 
j=i 

do201k=l / np 

j=j+l 

a(i, j)=x(i,l) **k 
201 continue 
n=np+l 
goto3 0 1 
c 
c 

c ************************************ 

400 continue 


1 go to 400 

*********************************** 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


if nd is not equal to 2 go to 500 
if (nd . ne . 2 ) goto500 

*********************************************************** 

*********************************************************** 

if we got here we have 2 design variables 

xl=x(i, 1) 
x2=x(i, 2) 

********************* 

add the constant and linear terms 

a(i,l)=l. 
a(i,2)=xl 
a (i , 3 ) =x2 
n=3 

if (np. It . 2) got o3 01 

******************** 

add the 2nd order terms 

a ( i , 4 ) =xl**2 
a ( i , 5) =xl*x2 
a ( i , 6) =x2**2 
n=6 

if (np. It . 3 ) goto301 


*★* 
* * * 
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c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 


******************* 


add the cubic terms 

a ( i , 7 ) =xl **3 
a(i, 8 )=xl** 2 *x 2 
a(i, 9) =xl*x 2**2 
a (i, 10 ) =x 2**3 
n=10 

if (np.lt. 4)goto301 
******************* 


add the 4th order terms 

a (i, ll)=xl**4 
a(i f 12)=xl**3*x2 
a (i, 13)=xl**2*x2**2 
a ( i , 14 ) =xl*x2**3 
a ( i , 15) =x2**4 
n=15 

if (np. It . 5) goto301 

******************* 

add the 5th order terms 

a (i , 16) =xl**5 
a(i, 17)=xl**4*x2 
a(i, 18 ) =xl**3 *x2 **2 
a(i, 19) =xl**2*x2**3 
a (i , 20 ) =xl*x 2 **4 
a (i , 21 ) =x 2**5 
n =21 

if (np. It. 6)goto301 
****************** 


algorithm 


not programed 


for polynomials of order 


larger than 5 


write( 6 , *) ' 
write ( 6 , *) ' 
write (7 , *) ' 
write ( 7 , *) ' 
stop 


for two design 
polynomials of 
for two design 
polynomials of 


variables, algorithm 
order larger than 5 ' 
variables, algorithm 
order larger than 5 ' 


not 

not 


programed for' 
programed for' 






** 
* * 


500 continue 


c if we got here number of design variables >2 

c ******************** 


enter constant and linear terms 

a ( i , 1 ) = 1 . 
j=l 
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do501)c=l , nd 
j=j+l 

a(i, j)=x(i,k) 

501 continue 
n=j 

if (np. It.2)goto301 

******************** 

enter the quadratic terms 

do502k=l,nd 
do502L=k, nd 

j=j+l 

a(i, j)=x(i,k)*x(i,L) 

502 continue 
n=j 

if (np. It . 3) goto301 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


******************** 


algorithm not programmed for more than quadratic approximation 
when number of design variables >2 


write (6 , *) ' 
write (6 , *) ' 
write (7,*) ' 
write (7 , *) ' 
stop 


ilgorithm not programmed for more than quadratic 
ipproximation when number of design variables >2 
ilgorithm not programmed for more than quadratic 
ipproximation when number of design variables >2 






*** 

**★ 


********** 

********* 


print out some results 

301 continue 

if (ip.lt. 4)goto302 

write (6,*) ' a(i,j)',(a(i,j) ,1-1, n) 
write (6/*) ' ' ...... 

write (7,*) ' a ( i , j ) ' , (a (l , D ) , 1=1 , n ) 
write(7 , *) ' 9 

302 continue 


c 

c 

c 


******** 




300 continue - , _ 

write(6,*)' number of undetermined coef-n- ,n 
write(7,*)' number of undetermined coef^-'^n 


c 


C 


return 

end 

SUBROUTINE PSI ( IP , M , N , DUMA , Y , XX) 
IMPLICIT REAL* 8 (A-H,0-Z) 


iMENSION DUMa(250, 136) r^x/oi oil RPTf21 21^ 

PENSION A ( 2 1 , 21 ) ,B(21,21) f D ( 2 i , 2 1) ,DI ( 21 , 2 1) , BPI ( 2 1 , 1 ) 

IMENSION C(21,21) , FI (21,21) , CPI (21,21) ,H (21,21) , HI (21, 21) 
IMENSION API (21, 21) 

IMENSION F ( 2 1 , 21 ) 

IMENSION IPIVOT (21), IWK (21,2) 
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DIMENSION y (250) 

DIMENSION XX (21) 

C 

c THIS SUBROUTINE CALCULATES PSEUDO INVERSE OF MATRIX A 

C M = ROW DIMENSION OF A LESS THAN N 

C N = COLUMN DIMENSION OF A 

C 

C COPY DUMA TO A 

C 

D090I=1 , M 
DO90J=l,N 

90 A (I , J) =DUMA ( I , J) 

C 

C 

C PRINT MATRIX A 

if (ip. It . 4) goto50 
WRITE (6, *) ' MATRIX A' 

WRITE (7,*)' MATRIX A' 

CALL WRITIT (M , N , A ) 

50 continue 
C 

C SET UP MATRIX B 

C 

D0100I=1,M 
D0100J=1,M 
100 B(I,J)=A(I,J) 

if (ip.lt. 4)goto51 
WRITE (6 , *) ' MATRIX B' 

WRITE (7,*)' MATRIX B' 

CALL WRITIT (M,M,B) 

51 continue 
C 

C GET D= B TRAN * B 

C 


DO101I=l,M 

doioij=i,m 

D (I , J) =0 . 

DO101K=l,M 

101 D (I , J) =D (I , J) +B (K, I) *B (K, J) 
if (ip. lt.4)goto52 
WRITE ( 6 , * ) ' MATRIX D' 

WRITE (7,*)' MATRIX D' 

CALL WRITIT (M,M,D) 

52 continue 

GET INVERSE OF D=DI 

MAX=21 

MDUM=0 

IOP=0 


300 


53 


CALL MATINV (MAX,M, D 
WRITE ( 6 , * ) ' DETERM= 
WRITE {!,*)' DETERM= 
D0300I=1 ,M 
D0300J=1 , M 


MDUM , DI , IOP , DETERM 
, DETERM, ' ISCALE=' 

, DETERM , ' ISCALE = ' 


ISCALE, IPIVOT, IWK) 

I SCALE 

ISCALE 


DI (I , J) =D (I , J) 
if (ip. It. 4) goto53 
WRITE ( 6 , * ) ' MATRIX DI' 
WRITE (7 , *) ' MATRIX DI ' 
CALL WRITIT (M, M, DI) 
continue 
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c GET PSEUDO INVERSE OF B = BPI = DI * B TRANS 

C 

DO102I=l ,M 
DO102 JQ=1,M 
BPI (I, JQ)=0. 

DO102 J=1 , M 

102 BPI (I, JQ)=BPI (I, JQ) +DI (I , J) *B(JQ, J) 
if ( ip . It . 4 ) goto54 

WRITE (6,*)' MATRIX BPI' 

WRITE (7,*)' MATRIX BPI' 

CALL WRITIT (M , M, BPI ) 

54 continue 
C 

C SET UP MATRIX C = A 

C 

DO103I=l , M 
D0103 J=1 , N 

103 C (I , J) =A (I , J) 

if (ip. lt.4)goto55 
WRITE (6,*)' MATRIX C' 

WRITE (7 , * ) ' MATRIX C' 

CALL WRITIT (M,N,C) 

55 continue 
C 

C SET UP MATRIX F = C * C TRANS 

C 

DO104I=l,M 
D0104 J=1 , M 
F(I,J)=0. 

D0104K=1,N 

104 F (I , J) =F ( I , J) +C (I,K) *C ( J , K) 
if (ip. It . 4) goto56 

WRITE ( 6 , * ) ' MATRIX F' 

WRITE (7,*)' MATRIX F' 

CALL WRITIT (M,M,F) 

56 continue 
C 

C GET THE INVERSE OF F = FI 

CALL MATINV (MAX , M , F , MDUM , FI , IOP , DETERM, I SCALE, IPIVOT, IWK) 
WRITE (6,*)' DETERM= ', DETERM , ' ISCALE=' , I SCALE 
WRITE(7,*)' DETERM= DETERM, ' ISCALE=' , ISCALE 
D03 011=1 , M 
D0301J=1,M 
301 FI (I , J) =F ( I , J) 

if (ip. lt.4)goto57 
WRITE ( 6 , * ) ' MATRIX FI' 

WRITE (7,*)' MATRIX FI' 

CALL WRITIT (M,M, FI) 

57 continue 

c 

c GET THE PSEUDO INVERSE OF C = CPI = C TRANS * FI 

C 

D0105IQ=1 , N 
D0105J=1,M 
CPI (IQ, J) =0 . 

D0105I=1,M 

105 CPI (IQ, J)=CPI(IQ, J)+C(I,IQ)*FI(I, J) 
if (ip. lt.4)goto58 
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WRITE (6 , * ) ' MATRIX CPI ' 

WRITE (7,*)' MATRIX CPI' 

CALL WRITIT (N , M, CPI ) 

58 continue 
C 

C SET UP MATRIX H = PSEUDO INVERSE OF B = BPI 


C 

C 


C 

C 

C 


DOl 061=1 , M 
D0106 J=1 , M 
106 H (I , J) =BPI ( I , J) 

if (ip.lt. 4)goto59 
WRITE (6,*)' MATRIX H' 
WRITE ( 7 , * ) ' MATRIX H' 
CALL WRITIT (M, M, H) 

59 continue 


GET INVERSE OF H = HI 
CALL MATINV (MAX , M , H , MDUM , HI , 
WRITE ( 6 , * ) ' DETERM=' , DETERM, 
WRITE (7,*)' DETERM=' , DETERM, 
D0302I=1 , M 
D0302 J=1 , M 
302 HI (I , J) =H (I , J) 

if (ip. It. 4) goto60 
WRITE (6,*)' MATRIX HI' 

WRITE ( 7 , * ) ' MATRIX HI' 

CALL WRITIT (M,M, HI) 

60 continue 


IOP, DETERM, ISCALE, IPIVOT, IWK) 
' I S CALE=' , ISCALE 
' ISCALE=' , ISCALE 


GET PSEUDO INVERSE OF A = API = CPI * HI * BPI 


DO107IQ=l , N 
D0107 J=1 , M 
API(IQ, J)=o. 

DO107I=l,M 

DO107K=l,M 

107 =API (IQ, J) ++CPI (IQ, I) *HI ( I , K) *BPI (K, J) 

if (ip.lt. 4)goto61 \ > 

WRITE (6, *) ' MATRIX API' 

WRITE (7,*)' MATRIX API' 

CALL WRITIT (N,M, API) 

61 continue 
C 

C GET XX = API * Y 
C 


D0108 IQ=1 , N 
XX (IQ) =0 . 

DO108 J=1 , M 

108 XX(IQ)=XX(IQ)+API(IQ,J)* Y (J) 
JDUM=1 

if (ip. It . 4) goto62 
WRITE (6,*)' MATRIX XX' 

WRITE (7 , *) ' MATRIX XX' 

CALL WRITIT (N,JDUM, XX) 

62 continue 


RETURN 

END 

SUBROUTINE WRITIT (MM, NN , XX) 
IMPLICIT REAL* 8 (A-H,0-Z) 
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DIMENSION XX (2 1,1) 

1 FORMAT (IX) 

2 FORMAT (10F7 . 2) 

WRITE (6, 1) 

DO100I=l , MM 

WRITE (6,2) (XX(I,J) ,J=1,NN) 

WRITE (7,2) (XX(I,J) ,J=1,NN) 

100 CONTINUE 
RETURN 
END 

SUBROUTINE EXACT ( IP , M, A, Y , B) 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION a (250, 136) ,b(136) ,y (250) 

DIMENSION IPIVOT(250) ,IWK(250,2) 

DIMENSION C (136 , 1) 

D0100I=1,M 

100 C (I , 1) =Y ( I ) 

MAX=250 

MDUM=1 

IOP=0 

CALL MATINV (MAX ,M, A, MDUM , C , IOP , DETERM , ISCALE , IPIVOT , IWK) 
WRITE (6,*)' DETERM=' , DETERM, ' ISCALE= ISCALE 
WRITE (7,*)' DETERM= ' , DETERM , ' I SCALE= ISCALE 
DO101I=l , M 
B ( I ) =C (1,1) 

101 CONTINUE 

if (ip. It . 4) goto50 

WRITE ( 6 , * ) ' MATRIX B' , (B ( I ) , 1=1 , M) 

WRITE (7 , *) ' MATRIX B' , (B ( I ) , 1=1 , M) 

50 continue 
RETURN 
END 

SUBROUTINE OVER ( IP , M, N, A, Y , B) 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION a (250, 136) ,b(136) ,y(250) 

DIMENSION IPIVOT (136) , IWK (136,2) 

DIMENSION AT A (136, 136) ,ATY(136,1) 

D0200I=1,N 
D0200J=1,N 
ATA(I , J) =0 . 

D0200K=1,M 

200 ATA(I, J) =ATA(I, J) +A(K, I) *A(K, J) 

DO201I=l , N 

ATY ( I , 1) =0 . 

DO201K=l,M 

201 ATY (1,1) =ATY (1,1) +A (K , I ) *Y (K) 

MAX= 136 

MDUM=1 

IOP=0 

CALL MATINV (MAX , N , ATA , MDUM , ATY , IOP , DETERM , ISCALE , IPIVOT , IWK ) 

WRITE (6,*)' DETERM=' , DETERM, ' ISCALE= ISCALE 

WRITE (7,*)' DETERM=' , DETERM, ' ISCALE=' , ISCALE 

D0101I=1,N 

B ( I ) =ATY (1,1) 

101 CONTINUE 

if (ip. It.4)goto50 

WRITE (6, *) ' MATRIX B' , (B(I) , 1=1, N) 

WRITE(7 ,*) ' MATRIX B ' , (B ( I ) , 1=1 , N) 

50 continue 
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o o o 


RETURN 

END 

!:pucii n ?e^r 

I sis sesse sissrss r«y y ,^ ™ es 

dimension y(250) ,yhat(250) 

dimension yy (2000) , yyhat (2000) 
yo — o • 

dol00id=l,m 
yb=yb+y (id) 

100 continue 

yb=yb/f loat (m) 
error= 0 . 
doloiid=l,m 

, ni err °^ =error + (y ( id) -yhat(id ) ) **2 

tloat (“) ) /yb* ( 100 . ) 

writer?'*!/ error over designs=error = '.error 
avera ? e Y over design = yb = ' yb 
write(6,*)' coefficient v (as %)= / i ,Y 
write ( 7 ,*), coefficient v fas %)= .\l 

dd=0. 

do7769id=l , m 
dn=dn+ (yhat ( id) -yb) **2 
dd=dd+(y (id) -yb) **2 
7769 continue 

r2=dn/dd* ( 100 . ) 

wr ite(6, *) ' coefficient r 2 (as%^ 

c "e r rj' 7 '*>' r2 S3} : .'M 

perror= 0 . 
yg=o. 

dol55id=l , ng 
yg=yg+yy(id) 

155 SnS err0r+(yy ( id) - yyhat <ia» 0*2 
yg=yg/fioat(ng) 

writeT? ( S ff 3 r ° r/f : 1 ° at (n?) ) /Y9* ( 100 . ) 

write (7^*) ' Iverage r y S ove? 1 grid U - ^ yq =P ^ ^r0r=, ' perror 

write (6 , *) • coefficient™ f < ~ Y9 ' yg 

write (7,*)' coefficient vg = /'va 
return y ' vg 

end 

implicit re^l 8 N yi^o- 2 ^' M,B,I ° P ' DETERM,ISCALE ' IPIVOT ^WK) MATINV 2 

PURPOSE - MATINV INVERTS A REAL SQUARE MATRIX A MATINV 5 

IN ADDITION THE ROUTINE SOLVES^THE MATRIX MATINV 6 

mAlRIX MATINV 7 
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c 

c 

c 

c 

c 

c 

c 

c USE 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


EQUATION AX=B f^I S^ALSC) 1 AN^OPTIO^T^HAVE^ THE^ 

?o L S^u^f f 

FOR SAVING TIME AND STORAGE. 

CALL MATIHV < MAX , N , A , M , B , IOP , DETERM , ISCALE , I PIVOT , IWK) 
MAX 


™e S ?S“p™oram . 


N 

A 


- the ORDER OF A, l.LE.N.LE.MAX. 


M 


at least n. 

the NUMBER OF COLUMN VECTORS IN ?c 
m=0 SIGNALS THAT THE SUBROUTINE IS 
USED SOLELY FOR INVERSION, HOWEVER, 

m ?h! ca£l statement an entry corre- 

SPONDING TO B MUST BE PRESENT. 

. A two-dimensional array of the constant 

VFfTOR B ON RETURN TO CALLING PROGRAM , 

X S STORED IN B. B SHOULD HAVE ITS FIRST 
dwInskn MAX AND ITS SECOND AT LEAST M. 

COMPUTE DETERMINANT OPTION. 

IOP=0 COMPUTES THE MATRIX INVERSE AND 
determinant . 

IOP=l COMPUTES THE MATRIX INVERSE ONLY. 

°de£ iwm) < ““Too ( ISCALE, , 

THE COMPUTATION ° ET ( A ) SHOULD SINCE IF 
ATTEMPTED IN THE USER PROGRAM SINCE IF 
THE ORDER OF A IS LARGER AND/OR THE 
MAGNITUDE OF ITS ELEMENTS ARE LARGE ( SMALL) , 
WE DE?( A) CALCULATION MAY CAUSE OVERFLOW 
(UNDERFLOW) ! DETERM SET TO ZERO FOR 
qTNGULAR MATRIX CONDITION, FOR EITHER 
TOP=1?5 r 0? should be checked BY PROGRAMER 
ON RETURN TO MAIN PROGRAM. 


B 


IOP - 


ISCALE - 


IPIVOT - 


A SCALE FACTOR COMPUTED BY THE 
SUBROUTINE TO AVOID OVERFLOW OR 
UNDERFLOW IN THE COMPUTATION OF 
THE QUANTITY , DETERM . 

A ONE DIMENSIONAL INTEGER ARRAY 
USED BY THE SUBPROGRAM TO STORE 
PIVOTOL INFORMATION. IT SHOULD BE 
DIMENSIONED AT LEAST N. IN GENERAL 


MATINV 8 
MATINV 9 
MATINV10 

matinvii 

MATINV 12 

MATINV13 

MATINV14 

MATINV15 

MATINV16 

MATINV 17 

MATINV 18 

MATINV19 

MATINV 2 0 

MATINV 2 1 

MATINV22 

MATINV 2 3 

MATINV24 

MATINV 2 5 

MATINV 2 6 

MATINV27 

MATINV28 

MATINV29 

MATINV30 

MATINV 3 1 

MATINV 3 2 

MATINV 3 3 

MATINV34 

MATINV 3 5 

MATINV 3 6 

MATINV37 

MATINV38 

MATINV39 

MATINV40 

MATINV4 1 

MATINV42 

MATINV4 3 

MATINV 4 4 

MATINV4 5 

MATINV46 

MATINV47 

MATINV48 

MATINV49 

MATINV50 

MATINV51 

MATINV52 

MATINV53 

MATINV54 

MATINV55 

MATINV56 

MATINV57 

MATINV58 

MATINV59 

MATINV60 

MATINV61 

MATINV62 

MATINV63 

MATINV64 

MATINV65 

MATINV66 

MATINV67 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THE USER DOES NOT NEED TO MAKE USE 
OF THIS ARRAY. 

IWK - A TWO-DIMENSIONAL INTEGER ARRAY OF 

TEMPORARY STORAGE USED BY THE ROUTINE. 
IWK SHOULD HAVE ITS FIRST DIMENSION 
MAX, AND ITS SECOND 2. 


REQUIRED ROUTINES- 
REFERENCE 


STORAGE 


-FOX,L, AN INTRODUCTION TO NUMERICAL 
LINEAR ALGEBRA 

- 542 OCTAL LOCATIONS 


LANGUAGE -FORTRAN 

LIBRARY FUNCTIONS -ABS 


RELEASED 
LATEST REVISION 


- JULY 1973 


- JULY 29, 1981 

COMPUTER SCIECES CORPORATION 
HAMPTON , VA 

C* *************************************************************** 

c 

DIMENSION IPIVOT(N) ,A(MAX,N) ,B(MAX,N) , IWK (MAX, 2) 

EQUIVALENCE (IROW,JROW), ( ICOLUM , JCOLUM) , (AMAX, T, SWAP) 


20 


C 

C 

C 


INITIALIZATION 

ISCALE=0 

Rl= ( 10 . 0d+00) **32 
R2=l . 0d+00/Rl 
DETERM=1 . 0d+00 
DO 20 J=1 , N 
IPIVOT ( J) =0 
CONTINUE 
DO 550 1=1, N 

SEARCH FOR PIVOT ELEMENT 

AMAX=0. 0d+00 
DO 105 J=1,N 

IF (IPIVOT ( J) -1) 60, 105, 60 


C 

c 


60 

DO 100 K=1,N 



IF ( IPIVOT (K) -1 ) 80, 

100, 

80 

TMAX = ABS (A ( J, K) ) 



IF (AMAX-TMAX) 85,100, 

r 100 

85 

IROW=J 

ICOLUM=K 

AMAX=TMAX 


100 

CONTINUE 


105 

CONTINUE 

IF (AMAX) 740,106,110 


106 

DETERM=0. 0d+00 
ISCALE=0 
GO TO 740 


110 

IPIVOT (ICOLUM) = 1 



INTERCHANGE ROWS TO PUT PIVOT 


MAT I NV 6 8 
MATINV69 
MATINV70 
MATINV7 1 
MATINV7 2 
MATINV73 
MATINV74 
MATINV75 
MATINV76 
MATINV77 
MATINV7 8 
MATINV7 9 
MATINV80 
MATINV8 1 
MATINV82 
MATINV83 
MATINV84 
MATINV85 
MATINV86 
MAT I NV 8 7 
MATINV88 
MATINV89 
MATINV90 
****** *MATINV9 1 
MATINV92 
MATINV93 
MATINV94 
MATINV98 
MATINV99 
MATIN100 
MATIN101 
MATIN102 
MATIN103 
MATIN104 
MATIN105 
MATIN106 
MATIN107 
MATIN108 
MATIN1 09 
MATIN1 10 
MATIN111 
MATIN1 12 
MATIN113 
MATIN1 14 
MATIN115 
MATIN1 16 
MATIN1 17 
MATIN118 
MATIN119 
MATIN12 0 
MATIN12 1 
MATIN12 2 
MATIN1 2 3 
MATIN12 4 
MATIN125 
MATIN12 6 
MATIN12 7 
MATIN128 
MATIN129 
MATIN1 3 0 
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c 

140 


200 

210 

250 

260 


C 

C 

C 

1000 

1010 


1020 


1030 

1040 


1050 

1060 

1070 


1080 


1090 

2000 


2010 

320 
C 

C 

C 

321 

350 

360 

370 

C 

C 


IF (IROW-ICOLUM) 140, 260, 140 

DETERM=-DETERM 
DO 200 L=1,N 
SWAP=A ( IROW , L) 

A ( IROW , L) =A ( ICOLUM , L) 

A ( ICOLUM , L) =SWAP 
CONTINUE 

IF (M) 260, 260, 210 
DO 250 L=l, M 
SWAP=B ( IROW , L) 

B ( IROW, L)=B( ICOLUM, L) 

B (ICOLUM, L)=SWAP 
CONTINUE 
IWK (1,1) =IROW 
IWK (1,2) =ICOLUM 
PIVOT=A( ICOLUM, ICOLUM) 

IF(IOP) 740,1000,321 

SCALE THE DETERMINANT 

PIVOTI=PIVOT 

IF (ABS( DETERM) -Rl) 1030,1010,1010 

DETERM=DETERM/R1 

ISCALE=I SCALE+ 1 

IF (ABS (DETERM) -Rl) 1060,1020,1020 
DETERM=DETERM/R1 
ISCALE=ISCALE+1 
GO TO 1060 

IF (ABS (DETERM) -R2 ) 1040,1040,1060 

DETERM=DETERM*R1 

ISCALE=ISCALE-1 

IF (ABS (DETERM) -R2) 1050,1050,1060 

DETERM=DETERM*R1 

ISCALE=ISCALE-1 

IF (ABS (PIVOTI ) -Rl) 1090,1070,1070 
PIVOTI=PIVOTI /Rl 
ISCALE=ISCALE+1 

IF (ABS ( PIVOTI )-Rl) 320,1080,1080 
PIVOTI=PIVOTI /Rl 
ISCALE=ISCALE+1 
GO TO 320 

IF (ABS (PIVOTI) -R2) 2000, 2000, 320 

PIVOTI=PIVOTI*Rl 

ISCALE=ISCALE-1 

IF (ABS (PIVOTI) -R2 ) 2010,2010,320 
PI VOTI=PI VOTI *R1 
ISCALE=ISCALE-1 
DETERM=DETERM* PIVOTI 

DIVIDE PIVOT ROW BY PIVOT ELEMENT 

A (ICOLUM, ICOLUM) =1 . 0d+00 
DO 350 L=1,N 

A (ICOLUM, L) =A ( ICOLUM, L) /PIVOT 
IF (M) 380, 380, 360 

DO 370 L=1 , M 

B( ICOLUM, L)=B( ICOLUM, L) /PIVOT 
REDUCE NON-PIVOT ROWS 


MATIN13 1 

MATIN132 

MATIN133 

MATIN134 

MATIN135 

MATIN13 6 

MATIN137 

MATIN138 

MATIN139 

MATIN140 

MATIN14 1 

MATIN14 2 

MATIN14 3 

MATIN144 

MATIN14 5 

MATIN14 6 

MATIN147 

MATIN148 

MATIN149 

MATIN150 

MATIN151 

MATIN152 

MATIN153 

MATIN154 

MATIN155 

MATIN156 

MATIN157 

MATIN158 

MATIN159 

MATIN160 

MATIN16 1 

MATIN162 

MATIN163 

MATIN164 

MATIN165 

MATIN166 

MATIN167 

MATIN168 

MATIN169 

MATIN17 0 

MATIN17 1 

MATIN172 

MATIN173 

MATIN 17 4 

MATIN17 5 

MATIN17 6 

MATIN177 

MATIN17 8 

MATIN17 9 

MATIN180 

MATIN181 

MATIN182 

MATIN18 3 

MATIN184 

MATIN185 

MATIN186 

MATIN187 

MATIN188 

MATIN189 

MATIN190 
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non 


c 

380 DO 550 L1=1,N 

IF(Ll-ICOLUM) 400 , 550, 400 
400 T=A( LI, ICOLUM) 

A (LI , ICOLUM) =0 . 0d+00 
DO 450 L=1,N 

450 A ( L 1 , L ) = A ( L 1 , L ) -A ( ICOLUM, L) *T 

IF (M) 550, 550, 460 

460 DO 500 L=1 , M 

500 B (LI , L) =B (LI , L) -B ( ICOLUM, L) *T 

550 CONTINUE 

INTERCHANGE COLUMNS 


C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DO 710 1=1, N 
L=N+1-I 

IF ( IWK (L, 1) -IWK(L, 2) ) 630, 710, 630 
630 JROW=IWK(L, 1) 

J COLUM=I WK ( L , 2 ) 

DO 705 K=1 , N 
SWAP=A(K, JROW) 

A (K , JROW) =A (K , JCOLUM) 

A (K , JCOLUM) =SWAP 
705 CONTINUE 
710 CONTINUE 
740 RETURN 
END 

ROUTINE NAME - HC318=EPSLON 
FROM EISPACK 


LATEST REVISION 

- AUGUST 1,1984 

COMPUTER SCIENCES CORP. , HAMPTON, VA. 

PURPOSE 


- THE FORTRAN FUNCTION EPSLON ESTIMATES UNIT 
ROUNDOFF IN QUANTITIES OF SIZE X. 

USAGE 


- VARIABLE = EPSLON (X) 

ARGUMENTS 

X 

- IS A REAL INPUT VARIABLE WHICH REPRESENTS THE 
QUANTITIES OF SIZE IN WHICH UNIT ROUNDOFF 
WILL BE ESTIMATED. 

REQUIRED ROUTINES 

- NONE 


REMARKS 1. IT SHOULD BE NOTED THAT EPSLON IS A FUNCTION 
DESIGNED TO BE CALLED BY ROUTINES IN THE 
EISPACK VERSION 3. 


THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL 
SYSTEMS SATISFYING THE FOLLOWING TWO 
ASSUMPTIONS, 

A. THE BASE USED IN REPRESENTING FLOATING 
POINT NUMBERS IS NOT A POWER OF THREE. 

B. THE QUANTITY A IN STATEMENT 10 IS 
REPRESENTED TO THE ACCURACY USED IN FLOATING 
POINT VARIABLES THAT ARE STORED IN MEMORY. 


MATIN19 1 
MATIN192 
MATIN193 
MATIN194 
MATIN195 
MATIN196 
MATIN197 
MATIN198 
MATIN199 
MATIN200 
MATIN2 01 
MATIN2 02 
MATIN2 03 
MATIN2 04 
MATIN2 05 
MATIN2 06 
MATIN2 07 
MATIN2 08 
MATIN2 09 
MATIN210 
MATIN211 
MATIN2 12 
MATIN2 1 3 
MATIN2 14 
MATIN2 15 
MATIN2 16 
MATIN2 17 
EPSLON 2 
EPSLON 3 
EPSLON 4 
EPSLON 5 
EPSLON 6 
EPSLON 7 
EPSLON 8 
EPSLON 9 
EPSLON1 0 
EPSLON11 
EPSLON12 
EPSLON13 
EPSLON14 
EPSLON15 
EPSLON16 
EPSLON17 
EPSLON18 
EPSLON19 
EPSLON2 0 
EPSLON2 1 
EPSLON2 2 
EPSLON2 3 
EPSLON24 
EPSLON25 
EPSLON2 6 
EPSLON27 
EPSLON28 
EPSLON2 9 
EPSLON3 0 
EPSLON3 1 
EPSLON3 2 
EPSLON3 3 
EPSLON34 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

C100 

c 

c 

c 

CA = 

c 


THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE 
INTENDED TO FORCE OPTIMIZING COMPILERS TO 
GENERATE CODE SATISFYING ASSUMPTION 2. 


UNDER THESE ASSUMPTIONS, 
THAT, 


IT SHOULD BE TRUE 


A IS NOT EXACTLY EQUAL TO FOUR-THIRDS, 

B HAS A ZERO FOR ITS LAST BIT OR DIGIT, 

C IS NOT EXACTLY EQUAL TO ONE, 

EPS MEASURES THE SEPARATION OF 1.0 FROM THE 
NEXT LARGER FLOATING POINT NUMBER. 


EXAMPLE : 

PROGRAM TR ( OUTPUT, TAPE6=OUTPUT) 
REAL X 
X = 4. 

A = EPSLON (X) 

WRITE ( 6 , 100) A 
FORMAT (5H0A = ,G22.14) 

STOP 
END 

OUTPUT : 

. 56843418860808 E- 13 




C*F45V1P0* 

REAL* 8 FUNCTION EPSLON 


(X) 


C 

REAL* 8 X 
REAL* 8 A , B , C , EPS 
A = 4.0E0/3.0E0 
10 B = A - 1 . 0E0 
C = B + B + B 
EPS = ABS ( C-l . 0E0 ) 

IF (EPS .EQ. 0.0E0) GO TO 10 
EPSLON = EPS*ABS (X) 

RETURN 

C** THIS PROGRAM VALID ON FTN4 AND FTN5 ** 


C 

C 

C 

C- 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 


END 

ROUTINE NAME 
FROM EISPACK 


LATEST REVISION 


PURPOSE 


- PF260=QZHES 


- AUGUST 1,1984 

COMPUTER SCIENCES CORP. , HAMPTON, VA. 


- THIS SUBROUTINE ACCEPTS A PAIR OF REAL 

GENERAL MATRICES AND REDUCES ONE OF THEM TO 
UPPER HESSENBERG FORM AND THE OTHER TO UPPER 
TRIANGULAR FORM USING ORTHOGONAL 
TRANSFORMATIONS. IT IS USUALLY FOLLOWED BY 
QZIT (PF261) , QZVAL(PF262) AND, POSSIBLY, 
QZVEC (PF2 63 ) . 


EPSLON35 
EPSLON36 
EPSLON 3 7 
EPSLON38 
EPSLON39 
EPSLON40 
EPSLON41 
EPSLON42 
EPSLON43 
EPSLON44 
EPSLON45 
EPSLON46 
EPSLON47 
EPSLON48 
EPSLON49 
EPSLON50 
EPSLON51 
EPSLON52 
EPSLON53 
EPSLON54 
EPSLON55 
EPSLON56 
EPSLON57 
EPSLON58 
EPSLON59 
EPSLON60 
EPSLON61 
EPSLON62 
EPSLON63 
EPSLON64 
EPSLON65 
EISPAK 
EISPAK32 
EISPAK 
EISPAK 
EISPAK35 
EISPAK36 
EISPAK37 
EISPAK38 
EISPAK39 
EISPAK4 0 
EISPAK41 
EISPAK42 
EISPAK43 
QZHES 2 
QZHES 3 
QZHES 4 
QZHES 5 
QZHES 6 
QZHES 7 
QZHES 8 
QZHES 9 
QZHES 10 
QZHES 11 
QZHES 12 
QZHES 13 
QZHES 14 
QZHES 15 
QZHES 16 
QZHES 17 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


USAGE 

ARGUMENTS 


NM 


N 


- CALL QZHES(NM,N,A^B,MATZ, Z) 

- ON INPUT NM MUST BE SET TO THE ROW DIMENSION 
OF TWO-DIMENSIONAL ARRAY PARAMETERS AS 
DECLARED IN THE CALLING PROGRAM DIMENSION 
STATEMENT. 

- ON INPUT N IS THE ORDER OF THE MATRICES. 

- ON INPUT A CONTAINS A REAL GENERAL MATRIX. 
MUST BE OF DIMENSION NM X N. 


18 

19 

20 
21 
22 
23 


QZHES 
QZHES 
QZHES 
QZHES 
QZHES 
QZHES 
QZHES 24 
QZHES 25 
QZHES 
QZHES 
QZHES 
QZHES 
QZHES 
QZHES 
QZHES 


MATZ 


ON OUTPUT A HAS BEEN REDUCED TO UPPER 
HESSENBERG FORM. THE ELEMENTS BELOW THE FIRSTQZHES 
SUBDIAGONAL HAVE BEEN SET TO ZERO. QZHES 

QZHES 

ON INPUT B CONTAINS A REAL GENERAL MATRIX. QZHES 
MUST BE OF DIMENSION NM X N. QZHES 

QZHES 

ON OUTPUT B HAS BEEN REDUCED TO UPPER QZHES 

TRIANGULAR FORM. THE ELEMENTS BELOW THE MAIN QZHES 
DIAGONAL HAVE BEEN SET TO ZERO. QZHES 

QZHES 

ON INPUT MATZ SHOULD BE SET TO .TRUE. IF THE QZHES 


RIGHT HAND TRANSFORMATIONS ARE TO BE 
ACCUMULATED FOR LATER USE IN COMPUTING 
EIGENVECTORS, AND TO .FALSE. OTHERWISE. 

ON OUTPUT Z CONTAINS THE PRODUCT OF THE RIGHT 
HAND TRANSFORMATIONS IF MATZ HAS BEfiN SET TO 
-TRUE. OTHERWISE, Z IS NOT REFERENCED. 

MUST BE OF DIMENSION NM X N. 


REQUIRED ROUTINES 
REMARKS 


- NONE 


1 . 


THIS SUBROUTINE IS THE FIRST STEP OF THE QZ 
ALGORITHM FOR SOLVING GENERALIZED MATRIX 
EIGENVALUE PROBLEMS, SIAM J. NUMER. ANAL. 10 
241-256(1973) BY MOLER AND STEWART. 


EXAMPLE : 

PROGRAM TQZHES (OUTPUT , TAPE6=OUTPUT) 
DIMENSION A ( 5 , 5 ) ,Z(5,5) , B ( 5 , 5 ) 
LOGICAL MATZ 

N = 5 
NM = 5 

MATZ = .TRUE. 


DATA A 


DATA B 


/10. ,2. ,3. ,2*1. ,2 
1 • § 1 • j 2 . , 1 . ,9. 


, 12 . , 1 . , 2 . , 1 . , 3 . , 1 , 

,3*1. ,-l. , 1. , 15. 


, 11 , 

/ 


/ 12 . , 1 . , 1. ,2. ,2*1. ,14. ,1. ,-l. ,1. ,-l. ,1. , 

16. ,-l. ,1. ,2. ,-l. ,-l. ,12. ,-l. ,3*1. , -1 . ,11, 


/ 


CALL QZHES (NM,N, A, B, MATZ, Z) 

WRITE (6, 100) ((A(I,J) ,1=1,5) ,J=1,5) , ( (B(I,J) ,1=1,5) ,J=1,5) , 


QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 


26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 
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c 

C100 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c — 


* ((Z(I,J) ,J=1,5) 
F0RMAT(1H ,5H A = /5(1H , 5 (G8 . 2 , 2X) / ) ) 

* 5H B = /5(1H , 5 (G8 . 2 , 2X) / ) 

* 5H Z = / 5 ( 1H , 5 (G8 . 2 , 2X) / ) ) 


STOP 

END 

OUTPUT 

A = 
-9.9 
-2.4 
.91 
-3.8 
2.7 
B = 
- 12 . 

2.3 

-.34 

-3.8 

2.5 


4.1 

11 . 

.26 

2.0 

-1.5 

0. 
16. 
-3.0 
. 80 
-1.4 


0. 

-3.0 

-13. 

1.7 

-.99 

0. 

0. 

- 12 . 

-1.5 

-1.5 


0. 

0. 

3 . 3 
- 11 . 

1.4 

0. 

0. 

0. 

- 10 . 

-1.5 


1.0 

0. 

0. 

.26 

0. 

.95 

0. 

-.14 

0 . 

. 87E-01 

- . 24E-01 

. 43 

0 . 

. 24E-01 

. 16 

.89 

0. 

-.96 

.26 

. 22E-01 


0. 

0. 

0 . 

2.6 

- 11 . 

0. 

0. 

0 . 

0 . 

-13. 

0. 

- . 70E— 01 
-.90 
.43 

- . 89E-01 


C 

C 


SUBROUTINE QZHES (NM, N , A, B , MATZ , Z) 

implicit real*8 (a-h,o-z) 

INTEGER I , J , K , L , N , LB , LI , NM , NK1 , NM1 , NM2 
REAL* 8 A (NM, N) , B (NM, N) , Z (NM, N) 

REAL* 8 R, S,T,U1,U2 ,V1,V2 , RHO 

LOGICAL MATZ 

IF (-NOT. MATZ) GO TO 10 

DO 3 J = 1, N 

DO 2 1 = 1, N 

Z(I,J) = 0.0E0 

2 CONTINUE 

Z(J,J) = 1.0E0 

3 C0NT: [ NUE . . REDUCE B TO UPPER TRIANGULAR FORM 
10 IF* (N -LE. 1) GO TO 170 

NM1 = N - 1 

DO 100 L = 1, NM1 
LI = L + 1 
S = 0.0E0 


DO 20 I = LI, N 

S = S + ABS ( B ( I , L) ) 

20 CONTINUE 

IF (S . EQ . 0.0E0) GO TO 100 

S = S + ABS (B (L, L) ) 


QZHES 78 
QZHES 79 
QZHES 80 
QZHES 81 
QZHES 82 
QZHES 83 
QZHES 84 
QZHES 85 
QZHES 86 
QZHES 87 
QZHES 88 
QZHES 89 
QZHES 90 
QZHES 91 
QZHES 92 
QZHES 93 
QZHES 94 
QZHES 95 
QZHES 96 
QZHES 97 
QZHES 98 
QZHES 99 
QZHES100 
QZHES101 
QZHES102 
QZHES103 
QZHES104 
QZHES105 
- QZHES106 
EISP6685 

EISP6686 
EISP6687 
EISP66 
EISP66 
EISP6690 
EISP6691 
EISP6692 
EISP6693 
EISP6694 
EISP6695 
EISP6696 
EISP6697 
EISP6698 
EISP6699 
EISP67 00 
EISP6701 
EISP67 02 
EISP6703 
EISP6704 
EISP67 05 
EISP6706 
EISP6707 
EISP67 08 
EISP6709 
EISP67 10 
EISP67 11 
EISP6712 
EISP67 13 
EISP6714 
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R = 0.0E0 


c 


DO 25 I = L, N 




B(I,L) = B ( I , L) / 
R = R + B ( I , L) **2 

S 


25 

CONTINUE 


c 


R = SIGN ( SQRT (R) ,B(L, 
B (L, L) = B (L, L) + R 
RHO = R * B (L, L) 

.L) ) 

c 





DO 50 J = LI, N 




T = 0.0E0 


c 


DO 30 I = L, N 




T = T + B ( I , L) 

* B (I , J) 


30 

CONTINUE 

c 


T = -T / RHO 


c 





DO 40 I = L, N 




B ( I , J) = B ( I , J) 

+ T * B ( I , L) 


40 

CONTINUE 

c 

50 

CONTINUE 


c 


DO 80 J = 1, N 




T = 0.0E0 


c 


DO 60 I = L, N 



60 

T = T + B ( I , L) 

* A(I,J) 


CONTINUE 

c 


T = -T / RHO 


c 





DO 70 I = L, N 



70 

A ( I , J) = A (I , J) 

+ T * B ( I , L) 


CONTINUE 

c 

80 

CONTINUE 


c 


B (L, L) = -S * R 


c 


DO 90 I = LI, N 




B ( I , L) = 0.0E0 



90 

CONTINUE 



c 

100 CONTINUE 

REDUCE A TO UPPER HESSENBERG FORM, WHILE 

KEEPING B TRIANGULAR 

IF (N .EQ. 2) GO TO 170 
NM2 = N - 2 
C 

DO 160 K = 1, NM2 
NK1 = NM1 - K 

C for L=N-1 STEP -1 UNTIL K+l DO — 

DO 150 LB = 1, NK1 
L = N - LB 
LI = L + 1 

c ZERO A (L+l , K) 


EISP6715 

EISP67 16 

EISP67 17 

EISP6718 

EISP6719 

EISP6720 

EISP672 1 

EISP6722 

EISP6723 

EISP6724 

EISP6725 

EISP6726 

EISP6727 

EISP6728 

EISP6729 

EISP6730 

EISP673 1 

EISP6732 

EISP6733 

EISP6734 

EISP673 5 

EISP673 6 

EISP6737 

EISP6738 

EISP6739 

EISP67 4 0 

EISP6741 

EISP6742 

EISP674 3 

EISP6744 

EISP6745 

EISP674 6 

EISP674 7 

EISP6748 

EISP6749 

EISP6750 

EISP6751 

EISP6752 

EISP6753 

EISP6754 

EISP6755 

EISP6756 

EISP67 57 

EISP67 58 

EISP6759 

EISP6760 

EISP6761 

EISP67 62 

EISP6763 

EISP6764 

EISP6765 

EISP6766 

EISP6767 

EISP6768 

EISP6769 

EISP6770 

EISP6771 

EISP6772 

EISP67 7 3 

EISP6774 
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s = ABS ( A ( L , K) ) + ABS (A ( LI , K) ) 

EISP6775 

EISP6776 


IF (S .EQ. 0.0E0) GO TO 150 

EISP6777 


U1 = A (L, K) / S 

EISP6778 


U2 = A ( LI , K) / S 

EISP6779 


R = SIGN ( SQRT (U1*U1+U2 *U2 ) ,U1) 

EISP6780 


VI = -(Ul + R) / R 

EISP67 8 1 


V2 = -U2 / R 

EISP6782 


U2 = V2 / VI 

EISP6783 

c 

DO 110 J = K, N 

EISP6784 


EISP6785 


T = A (L, J) + U2 * A ( LI , J) 

EISP6786 


A (L, J) = A (L, J) + T * VI 

EISP6787 


A (LI , J) = A ( LI , J) + T * V2 

EISP6788 

110 

CONTINUE 

EISP6789 

c 


EISP6790 


A (LI , K) = 0.0E0 

EISP6791 

c 

DO 120 J = L, N 

EISP6792 


EISP6793 


T = B (L, J) + U2 * B (LI , J) 

EISP6794 


B (L, J) = B (L, J) + T * VI 

EISP6795 


B (LI , J) = B (LI , J) + T * V2 

EISP6796 

120 

CONTINUE 

EISP6797 

c 

ZERO B (L+l , L) 

EISP6798 


S = ABS (B (LI , LI) ) + ABS (B (LI , L) ) 

EISP6799 


IF (S .EQ. 0.0E0) GO TO 150 

EISP6800 


Ul = B (LI , LI) / S 

EISP6801 


U2 = B(L1,L) / S 

EISP6802 


R = SIGN ( SQRT (U1*U1+U2 *U2 ) ,U1) 

EISP6803 


VI = -(Ul + R) / R 

EISP6804 


V2 = -U2 / R 

EISP68 05 


U2 = V2 / VI 

EISP6806 

c 

DO 130 1=1, LI 

EISP6807 


EISP6808 


T = B ( I , LI ) + U2 * B ( I , L) 

EISP6809 


B (I , LI) = B ( I , LI) + T * VI 

EISP68 10 


B ( I , L) = B ( I , L) + T * V2 

EISP68 11 

130 

CONTINUE 

EISP68 12 

C 


EISP6813 


B ( LI , L) = 0.0E0 

EISP6814 

C 

DO 140 I = 1, N 

EISP6815 


EISP68 16 


T = A ( I , LI ) + U2 * A ( I , L) 

EISP68 17 


A(I , LI) = A ( I , LI ) + T * VI 

EISP68 18 


A ( I , L) = A ( I , L) + T * V2 

EISP6819 

140 

CONTINUE 

EISP6820 

c 

IF (.NOT. MATZ) GO TO 150 

EISP682 1 


EISP682 2 

C 

DO 145 I = 1, N 

EISP6823 


EISP6824 


T = Z ( I , LI) + U2 * Z ( I , L) 

EISP682 5 


Z ( I , LI) = Z (I , LI) + T * VI 

EISP6826 


Z ( I , L) = Z ( I , L) + T * V2 

EISP6827 

145 

CONTINUE 

EISP6828 

C 


EISP6829 

150 

CONTINUE 

EISP6830 

C 


EISP6831 

160 

CONTINUE 

EISP6832 

C 


EISP683 3 

170 

RETURN 

EISP6834 

c** 

THIS PROGRAM VALID ON FTN4 AND FTN5 
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END 


C 

C 

C 

C- 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ROUTINE NAME 
FROM EISPACK 


- PF2 6 1=QZIT 


LATEST REVISION 


PURPOSE 


USAGE 

ARGUMENTS NM 

N 

A 


B 


EPS1 


AUGUST 1,1984 

COMPUTER SCIENCES CORP. , HAMPTON, VA , 


- THIS SUBROUTINE ACCEPTS A PAIR OF REAL 
MATRICES, ONE OF THEM IN UPPER HESSENBERG 
FORM AND THE OTHER IN UPPER TRIANGULAR FORM. 
IT REDUCES THE HESSENBERG MATRIX TO 

QUASI -TRIANGULAR FORM USING ORTHOGONAL 
TRANSFORMATIONS WHILE MAINTAINING THE 
TRIANGULAR FORM OF THE OTHER MATRIX. IT IS 
USUALLY PRECEDED QZHES(PF260) AND FOLLOWED 
BY QZVAL ( PF2 62 ) AND, POSSIBLY, QZVEC (PF2 63 ) . 

- CALL QZIT (NM,N, A, B, EPS1,MATZ, Z, I ERR) 

- ON INPUT NM MUST BE SET TO THE ROW DIMENSION 
OF TWO-DIMENSIONAL ARRAY PARAMETERS AS 
DECLARED IN THE CALLING PROGRAM DIMENSION 
STATEMENT. 

- ON INPUT N IS THE ORDER OF THE MATRICES. 

- ON INPUT A CONTAINS A REAL UPPER HESSENBERG 
MATRIX. 

MUST BE OF DIMENSION NM X N. 

ON OUTPUT A HAS BEEN REDUCED TO 
QUASI -TRIANGULAR FORM. THE ELEMENTS BELOW TH 
FIRST SUBDIAGONAL ARE STILL ZERO AND NO TWO 
CONSECUTIVE SUBDIAGONAL ELEMENTS ARE NONZERO. 


ON INPUT B CONTAINS A REAL UPPER TRIANGULAR 
MATRIX. 

MUST BE OF DIMENSION NM X N. 

ON OUTPUT B IS STILL IN UPPER TRIANGULAR 
FORM, ALTHOUGH ITS ELEMENTS HAVE BEEN ALTER! 
THE LOCATION B(N,1) IS USED TO STORE EPS1 
TIMES THE NORM OF B FOR LATER USE BY QZVAL 
QZVAL ( PF2 62 ) AND QZVEC (PF2 63 ) . 


* * '-*-*■'*'* WUUU X W Ub i L ru'l -L XL, ^ Zj J. 1 

NEGLIGIBLE ELEMENTS. EPS1 = 0.0 (OR NEGATIVE) QZIT 
MAY BE INPUT, IN WHICH CASE AN ELEMENT WILL BEQZIT 
NEGLECTED ONLY IF IT IS LESS THAN ROUNDOFF QZIT 
ERROR TIMES THE NORM OF ITS MATRIX. IF THE QZIT 
INPUT EPS1 IS POSITIVE, THEN AN ELEMENT WILL QZIT 
BE CONSIDERED NEGLIGIBLE IF IT IS LESS THAN QZIT 
EPS1 TIMES THE NORM OF ITS MATRIX. A POSITIVEQZIT 
VALUE OF EPS1 MAY RESULT IN FASTER EXECUTION, QZIT 
BUT LESS ACCURATE RESULTS. q ZIX 


EISP68 3 5 

QZIT 

2 

QZIT 

3 

QZIT 

4 

■ QZIT 

5 

QZIT 

6 

QZIT 

7 

QZIT 

8 

QZIT 

9 

QZIT 

10 

QZIT 

11 

QZIT 

12 

QZIT 

13 

QZIT 

14 

QZIT 

15 

QZIT 

16 

QZIT 

17 

QZIT 

18 

QZIT 

19 

QZIT 

20 

QZIT 

21 

QZIT 

22 

QZIT 

23 

QZIT 

24 

QZIT 

25 

QZIT 

26 

QZIT 

27 

QZIT 

28 

QZIT 

29 

QZIT 

30 

QZIT 

31 

QZIT 

32 

QZIT 

33 

QZIT 

34 

QZIT 

35 

2QZIT 

36 

QZIT 

37 

QZIT 

38 

QZIT 

39 

QZIT 

40 

QZIT 

41 

QZIT 

42 

QZIT 

43 

QZIT 

44 

QZIT 

45 

QZIT 

46 

QZIT 

47 

QZIT 

48 

QZIT 

49 

QZIT 

50 


51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
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o n o 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


MATZ 


I ERR - 


REQUIRED ROUTINES 
REMARKS 1 


- HC3 18=EPSLON 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C99 

C 

C 

C100 

c 

c 

c 

c 

c 


EXAMPLE : 

PROGRAM TQZIT( OUTPUT, TAPE6=OUTPUT) 
DIMENSION A(5,5) , B ( 5 , 5) ,2(5,5) 
LOGICAL MATZ 


N = 5 
NM = 5 
MATZ = 
EPS1 = 

DATA A 


DATA B 


. TRUE . 

0. 0E0 

/ 10 . ,2 

1. ,-l. 

/ 12 . ,1 

16 . , -1 


,3 

1 . 


QZIT 

ON INPUT MATZ SHOULD BE SET TO .TRUE. IF THE QZIT 
RIGHT HAND TRANSFORMATIONS ARE TO BE QZIT 

ACCUMULATED FOR LATER USE IN COMPUTING QZIT 

EIGENVECTORS, AND TO .FALSE. OTHERWISE. 

QZ 1 i 

• ON INPUT Z CONTAINS, IF MATZ HAS BEEN SET TO QZIT 
TRUE THE TRANSFORMATION MATRIX PRODUCED IN QZIT 
THE REDUCTION BY QZHES (PF260) , IF PERFORMED , QZIT 
OR ELSE THE IDENTITY MATRIX. IF MATZ HAS BEENQZIT 
SET TO .FALSE., Z IS NOT REFERENCED. QZIT 

MUST BE OF DIMENSION NM X N. 

ON OUTPUT Z CONTAINS THE PRODUCT OF THE 

RIGHT HAND TRANSFORMATIONS (FOR BOTH STEPS) IFQZIT 

MATZ HAS BEEN SET TO .TRUE.. QZIT 

ON OUTPUT I ERR IS SET TO n7TT 

ZERO FOR NORMAL RETURN. r ,„ mm itti 

J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED QZIT 
WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. 

QZl 1 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 


THIS SUBROUTINE IS THE SECOND STEP OF THE QZ 
ALGORITHM FOR SOLVING GENERALIZED MATRIX 
EIGENVALUE PROBLEMS , SIAM J. NUMER. ANAL. 10, 
241-256(1973) BY MOLER AND STEWART, AS 
MODIFIED IN TECHNICAL NOTE NASA TN 
D-7 3 05 (197 3) BY WARD. 


, , 2 * 1 , 

, 2 . , 1 


, 12 . 
r 3*1 


, 1 . , 2 
, , - 1 . 


,1.,3.,1.,11< 
1 . , 15 . / 


,-l.,2.,2*l. ,14. , 1 . , 1 m 1* /"I • ' 1 • ' 

,1. ,2. ,12. ,-l.,3*l., -1.,11- 


CALL QZHES ( NM , N , A , B , MATZ , Z ) 

CALL QZIT (NM, N , A , B , EPS1 , MATZ , Z , I ERR) 

WRITE (6,99) I ERR 

FORMAT ( 1H1 , 8H IERR = ,14) T 

WRITE ( 6 , 100) ( ( A ( I , J) ,1=1/5) , J=l, 5) , ((B(I,J) , 1-1/ 5 ) , , ) 

* ( ( Z ( I , J) ,1=1/5) , J=l,5) 

FORMAT (1H , 5H A = / 5 ( 1H , 5 (G8 . 2 , 2X) / ) ) 

* 5H B = / 5 ( 1H , 5(G8.2,2X) /) 

* 5H Z = / 5 ( 1H , 5 (G8 . 2 , 2X) / ) ) 

STOP 

END 


61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
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c 

OUTPUT : 




c 





c 

I ERR = 

0 



c 

A = 




c 

-15. 

-1.3 

0 . 

0 . 

c 

1.1 

7.4 

0 . 

0 . 

c 

1.5 

-1.5 

-16. 

0 . 

c 

-2.2 

.96 

1.0 

-10. 

c 

-2.6 

-.31 

1.2 

1.7 

c 

B = 




c 

-9.9 

0 . 

0 . 

0 . 

c 

-.29 

17. 

0 . 

0 . 

c 

1.3 

-2 . 1 

-14 . 

0 . 

c 

-1.9 

1.7 

.96 

-11. 

c 

-2.6 

-.32 

1.3 

2.1 

c 

Z = 




c 

.28 

- . 7 IE-01 

. 16 

-.24 

c 

.52 

-.24 

-.66 

. 48 

c 

.49 

. 56 

.49 

. 45 

c 

-.60 

. 48 

-.29 

.44 

c 

-.25 

-.63 

.45 

.57 


c 

c 

SUBROUTINE QZIT (NM, N, A, B, EPS1 , MATZ 
C 



QZIT 

121 


QZIT 

122 


QZIT 

123 


QZIT 

124 

0 . 

QZIT 

125 

0 . 

QZIT 

126 

0 . 

QZIT 

127 

0 . 

QZIT 

128 

-8.6 

QZIT 

129 


QZIT 

130 

. 31E-12 

QZIT 

131 

0 . 

QZIT 

132 

0 . 

QZIT 

133 

0 . 

QZIT 

134 

-13. 

QZIT 

135 


QZIT 

136 

-.91 

QZIT 

137 

64E-01 

QZIT 

138 

. 7 5E-01 

QZIT 

139 

-.38 

QZIT 

140 

- . 94E-01 

QZIT 

141 


QZIT 

142 


07TT 

143 

5836 

SRR) 

eispc 


implicit real*8 (a-h,o-z) 

INTEGER I, J,K,L,N, EN, K1 , K2 , LD, LL, LI , NA, NM, ISH, ITN, ITS,KM1,LM1, 
X ENM2 , I ERR , LORI , ENORN 

REAL*8 A (NM, N) , B (NM, N) , Z (NM , N) 

REAL* 8 R , S , T , A1 , A2 , A3 , EP , SH ,U1,U2,U3 ,V1,V2, V3 , AN I /All, 

X A12,A21,A22,A33,A34 / A43 / A44 # BNI,Bll f B12 # B22 # B33 f B34 , 

X B44 , EPSA, EPSB , EPS1 , ANORM, BNORM, EPSLON 

LOGICAL MATZ , NOTLAS 
I ERR = 0 

C COMPUTE EPSA, EPSB 

ANORM = 0.0E0 
BNORM = 0.0E0 
C 

DO 30 I = 1, N 
AN I = 0.0E0 

IF (I .NE. 1) ANI = ABS (A(I, 1-1) ) 

BN I = 0.0E0 
C 

DO 20 J = I, N 

ANI = ANI + ABS ( A ( I , J) ) 

BNI = BN I + ABS ( B ( I , J) ) 

20 CONTINUE 
C 

IF (ANI . GT . ANORM) ANORM = ANI 
IF (BNI .GT. BNORM) BNORM = BNI 
30 CONTINUE 
C 

IF (ANORM .EQ. 0.0E0) ANORM = 1.0E0 
IF (BNORM .EQ. 0.0E0) BNORM = 1.0E0 
EP = EPS1 

IF (EP .GT. 0.0E0) GO TO 50 


C USE ROUNDOFF LEVEL IF EPS1 IS ZERO 

EP = EPSLON ( 1. 0E0) 

50 EPSA = EP * ANORM 
EPSB = EP * BNORM 


EISP68 37 
EISP68 3 8 
EISP68 3 9 
EISP68 
EISP68 
EISP68 4 2 
EISP684 3 
EISP6844 
EISP684 5 
EISP684 6 
EISP6847 
EISP684 8 
EISP684 9 
EISP6850 
EISP68 5 1 
EISP6852 
EISP6853 
EISP6854 
EISP6855 
EISP6856 
EISP6857 
EISP68 58 
EISP6859 
EISP68 60 
EISP6861 
EISP6862 
EISP68 63 
EISP68 64 
EISP6865 
EISP68 66 
EISP68 67 
EISP6868 
EISP6869 
EISP687 0 
EISP687 1 
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on o no 


REDUCE A TO QUASI -TRIANGULAR FORM, WHILE EISP6872 

KEEPING B TRIANGULAR EISP6873 

LORI = 1 EISP6874 

ENORN = N EISP6875 

EN = N EISP6876 

ITN = 30*N EISP6877 

BEGIN QZ STEP EISP6878 

60 IF (EN .LE. 2) GO TO 1001 EISP6879 

IF (.NOT. MATZ) ENORN = EN EISP6880 

ITS = 0 EISP6881 

NA = EN - 1 EISP6882 

ENM2 = NA - 1 EISP6883 

70 ISH = 2 EISP6884 

CHECK FOR CONVERGENCE OR REDUCIBILITY . EISP6885 

FOR L=EN STEP -1 UNTIL 1 DO — EISP6886 

DO 80 LL = 1, EN EISP6887 

LM1 = EN - LL EISP6888 

L = LM1 + 1 EISP6889 

IF (L .EQ. 1) GO TO 95 EISP6890 

IF (ABS (A (L, LM1) ) .LE. EPSA) GO TO 90 EISP6891 

80 CONTINUE EISP6892 

C EISP6893 

90 A (L, LM1) = 0.0E0 EISP6894 

IF (L .LT. NA) GO TO 95 EISP6895 

C 1-BY-l OR 2-BY-2 BLOCK ISOLATED EISP6896 

EN = LM1 EISP6897 

GO TO 60 EISP6898 

C CHECK FOR SMALL TOP OF B EISP6899 

95 LD = L EISP6900 

100 LI = L + 1 EISP6901 

Bll = B (L, L) EISP6902 

IF (ABS (Bll) .GT. EPSB) GO TO 120 EISP6903 

B (L, L) = 0.0E0 EISP6904 

S = ABS (A (L, L) ) + ABS (A (LI , L) ) EISP6905 

U1 = A (L, L) / S EISP6906 

U2 = A (LI , L) / S EISP6907 

R = SIGN ( SQRT (U1*U1+U2 *U2 ) ,U1) EISP6908 

VI = -(Ul + R) / R EISP6909 

V2 = -U2 / R EISP6910 

U2 = V2 / VI EISP6911 

C EISP6912 

DO 110 J = L, ENORN EISP6913 

T = A (L, J) + U2 * A (LI , J) EISP6914 

A (L, J) = A (L, J) + T * VI EISP6915 

A (LI , J) = A (LI , J) + T * V2 EISP6916 

T = B (L, J) + U2 * B (LI , J) EISP6917 

B (L, J) = B (L, J) + T * VI EISP6918 

B (LI , J) = B (LI , J) + T * V2 EISP6919 

110 CONTINUE EISP6920 

C EISP692 1 

IF (L .NE. 1) A(L, LM1) = -A (L, LM1) EISP6922 

LM1 = L EISP692 3 

L = LI EISP6924 

GO TO 90 EISP6925 

120 All = A (L, L) / Bll EISP692 6 

A21 = A (LI , L) / Bll EISP6927 

IF (ISH .EQ. 1) GO TO 140 EISP6928 

C ITERATION STRATEGY EISP6929 

IF (ITN .EQ. 0) GO TO 1000 EISP6930 

IF (ITS .EQ. 10) GO TO 155 EISP6931 
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ooo 


C DETERMINE TYPE OF SHIFT 

B22 = B (LI , LI) 

IF ( ABS ( B2 2 ) .LT. EPSB) B22 = EPSB 
B33 = B (NA , NA) 

IF (ABS ( B3 3 ) . LT . EPSB) B33 = EPSB 

B44 = B(EN,EN) 

IF (ABS ( B44 ) .LT. EPSB) B44 = EPSB 

A33 = A (NA, NA) / B33 

A34 = A (NA, EN) / B44 

A43 = A ( EN , NA) / B33 

A44 = A ( EN , EN) / B44 

B34 = B (NA , EN) / B44 

T = 0.5E0 * ( A43 * B34 - A33 - A44) 

R = T * T + A3 4 * A43 - A33 * A44 
IF (R .LT. 0.0E0) GO TO 150 

C DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A 

ISH = 1 
R = SQRT(R) 

SH = -T + R 


S = -T - R 

IF (ABS(S-A44) .LT. ABS(SH-A44)) SH = S 

LOOK FOR TWO CONSECUTIVE SMALL 

SUB-DIAGONAL ELEMENTS OF A. 

FOR L=EN-2 STEP -1 UNTIL LD DO — 

DO 130 LL = LD, ENM2 
L = ENM2 + LD - LL 
IF (L .EQ. LD) GO TO 140 
LM1 = L - 1 
LI = L + 1 
T = A (L, L) 

IF (ABS (B (L, L) ) .GT. EPSB) T = T - SH * B (L, L) 

IF (ABS ( A (L, LM1 ) ) . LE . ABS (T/A (LI , L) ) * EPSA) GO TO 100 
130 CONTINUE 


C 

140 A1 = All - SH 
A2 = A21 

IF (L .NE. LD) A (L, LM1) = -A (L , LM1 ) 

GO TO 160 

C DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A 

150 A12 = A(L,L1) / B22 

A22 = A (LI , LI ) / B22 

B12 = B(L,L1) / B22 

A1 = ( (A3 3 - All) * ( A44 - All) - A34 * A43 + A43 * B34 * All) 

X / A21 + A12 - All * B12 

A2 = (A22 - All) - A21 * B12 - (A33 - All) - (A44 - All) 

X + A43 * B34 

A3 = A(L1+1,L1) / B22 
GO TO 160 

C AD HOC SHIFT 

155 A1 = 0.0E0 
A2 = 1.0E0 
A3 = 1.1605E0 
160 ITS = ITS + 1 
ITN = ITN - 1 
IF (.NOT. MATZ) LORI = LD 

C MAIN LOOP 

DO 260 K = L, NA 


NOTLAS = K .NE. NA .AND. ISH . EQ . 2 
K1 = K + 1 
K2 = K + 2 


EISP6932 

EISP693 3 

EISP6934 

EISP693 5 

EISP693 6 

EISP6937 

EISP6938 

EISP6939 

EISP6940 

EISP694 1 

EISP69 4 2 

EISP6943 

EISP69 44 

EISP6945 

EISP6946 

EISP694 7 

EISP6948 

EISP6949 

EISP6950 

EISP695 1 

EISP6952 

EISP6953 

EISP69 54 

EISP6955 

EISP6956 

EISP6957 

EISP6958 

EISP69 59 

EISP69 60 

EISP6961 

EISP6962 

EISP69 6 3 

EISP69 64 

EISP69 65 

EISP69 66 

EISP6967 

EISP6968 

EISP69 69 

EISP697 0 

EISP697 1 

EISP69 7 2 

EISP697 3 

EISP6974 

EISP697 5 

EISP697 6 

EISP6977 

EISP6978 

EISP69 7 9 

EISP6980 

EISP6981 

EISP6982 

EISP6983 

EISP6984 

EISP6985 

EISP6986 

EISP6987 

EISP6988 

EISP698 9 

EISP699 0 

EISP699 1 
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c 


170 


C 


180 

C 


C 

190 


200 


C 


210 

C 


C 

220 


KM1 = MAXO (K-l , L) 

LL = MINO (EN, Kl+ISH) 

IF (NOTLAS) GO TO 190 

ZERO A (K+l , K-l ) 

IF (K .EQ. L) GO TO 170 
A1 = A(K,KM1) 

A2 = A (K1 , KM1) 

S = ABS(Al) + ABS ( A2 ) 

IF (S .EQ. O.OEO) GO TO 70 
U1 = A1 / S 
U2 = A2 / S 

R = SIGN(SQRT(U1*U1+U2*U2) ,U1) 
VI = -(Ul + R) / R 
V2 = -U2 / R 
U2 = V2 / VI 

DO 180 J = KM1 , ENORN 

T = A (K, J) + U2 * A (K1 , J) 

A (K, J) = A (K, J) + T * VI 
A (K1 , J) = A (K1 , J) + T * V2 
T = B (K , J) + U2 * B (K1 , J) 

B (K, J) = B (K , J) + T * VI 
B (K1 , J) = B (K1 , J) + T * V2 
CONTINUE 


IF (K .NE. L) A (K1 , KM1 ) = O.OEO 
GO TO 240 

ZERO A (K+l , K-l ) AND A(K+2,K-1) 

IF (K .EQ. L) GO TO 200 
A1 = A (K, KM1 ) 

A2 = A (K1 , KM1 ) 

A3 = A (K2 , KM1) 

S = ABS(Al) + ABS ( A2 ) + ABS (A3) 

IF (S .EQ. O.OEO) GO TO 260 
Ul = A1 / S 

U2 = A2 / S 

U3 = A3 / S 

R = SIGN ( SQRT (U1*U1+U2 *U2+U3*U3 ) , Ul) 


VI = -(Ul + R) 
V2 = -U2 / R 
V3 = -U3 / R 
U2 = V2 / VI 
U3 = V3 / VI 


/ R 


DO 210 J = KM1 , ENORN 

T = A (K , J) + U2 * A (K1 , J) + U3 * A(K2,J) 
A (K , J) = A (K , J) + T * VI 
A (K1 , J) = A (K1 , J) + T * V2 

A (K2 , J) = A (K2 , J) + T * V3 

T = B (K , J) + U2 * B (K1 , J) + U3 * B(K2,J) 
B(K, J) = B (K , J) + T * VI 
B (K1 , J) = B (K1 , J) + T * V2 

B (K2 , J) = B (K2 , J) + T * V3 

CONTINUE 


IF (K .EQ. L) GO TO 220 
A (K1 , KM1) = O.OEO 
A (K2 , KM1 ) = O.OEO 

ZERO B (K+2 , K+l ) AND B(K+2,K) 

S = ABS (B (K2 , K2 ) ) + ABS (B (K2 ,K1) ) + ABS(B(K2,K)) 


EISP6992 

EISP6993 

EISP6994 

EISP6995 

EISP6996 

EISP6997 

EISP6998 

EISP6999 

EISP7000 

EISP7001 

EISP7002 

EISP7003 

EISP7004 

EISP7005 

EISP7006 

EISP7007 

EISP7008 

EISP7009 

EISP7010 

EISP7011 

EISP7012 

EISP7013 

EISP7014 

EISP7015 

EISP7016 

EISP7017 

EISP7018 

EISP7019 

EISP7020 

EISP7021 

EISP7022 

EISP7023 

EISP7024 

EISP7025 

EISP7026 

EISP7027 

EISP7028 

EISP7029 

EISP7030 

EISP703 1 

EISP703 2 

EISP703 3 

EISP7034 

EISP7035 

EISP7036 

EISP7037 

EISP7038 

EISP7039 

EISP7040 

EISP7041 

EISP704 2 

EISP704 3 

EISP7044 

EISP7045 

EISP704 6 

EISP7047 

EISP7048 

EISP7049 

EISP7050 

EISP7051 
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c 


230 

C 


C 


235 

C 

240 


C 


250 

C 


C 


255 

C 


IF (S .EQ. 0.0E0) GO TO 240 
U1 = B(K2,K2) / S 
U2 = B(K2,K1) / S 
U3 = B (K2 , K) / S 

R = SIGN (SQRT(U1*U1+U2*U2+U3*U3) , Ul) 


VI 

— 

" (ui 

+ R) / 

R 

V2 

= 

-U2 

/ R 


V3 

= 

-U3 

/ R 


U2 

= 

V2 / 

VI 


U3 

= 

V3 / 

VI 


DO 

230 I 

= LORI, 

LL 


T = A ( I , K2 ) + U2 * A ( I , K1 ) + U3 * A(I,K) 
A (I , K2 ) = A (I , K2 ) + T * VI 

A ( I , Kl) = A (I , Kl) + T * V2 

A ( I , K) = A ( I , K) + T * V3 

T = B (I , K2 ) + U2 * B (I , Kl) + U3 * B(I,K) 
B(I,K2) = B (I , K2 ) + T * VI 

B (I , Kl) = B (I , Kl) + T * V2 

B (I , K) = B ( I , K) + T * V3 
CONTINUE 


B (K2 , K) = 0.0E0 
B (K2 , Kl) = O.OEO 
IF (.NOT. MATZ) GO TO 240 

DO 235 1=1, N 

T = Z ( I , K2 ) + U2 * Z ( I , Kl) + U3 * Z(I,K) 
Z(I,K2) = Z ( I , K2 ) + T * VI 
Z(I,K1) = Z (I , Kl) + T * V2 
Z ( I / K) = Z ( I , K) + T * V3 
CONTINUE 

ZERO B (K+l , K) 

S = ABS (B (Kl , Kl) ) + ABS (B (Kl , K) ) 

IF (S .EQ. O.OEO) GO TO 260 
Ul = B (Kl , Kl ) / S 
U2 = B (Kl , K) / S 
R = SIGN (SQRT (U1*U1+U2 *U2 ) ,U1) 

VI = -(Ul + R) / R 
V2 = -U2 / R 
U2 = V2 / VI 

DO 250 I = LORI, LL 

T = A ( I , Kl ) + U2 * A ( I , K) 

A ( I , Kl ) = A ( I , Kl ) + T * VI 
A (I , K) = A ( I , K) + T * V2 
T = B ( I , Kl ) + U2 * B ( I , K) 

B ( I , Kl ) = B (I , Kl) + T * VI 
B ( I , K) = B ( I , K) + T * V2 
CONTINUE 


B (Kl , K) = O.OEO 
IF (.NOT. MATZ) GO TO 260 

DO 255 I = 1, N 

T = Z (I , Kl) + U2 * Z ( I , K) 

Z (I , Kl ) = Z ( I , Kl ) + T * VI 
Z ( I , K) = Z ( I , K) + T * V2 
CONTINUE 
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EISP7084 

EISP7085 

EISP7086 

EISP7087 

EISP7088 

EISP7089 

EISP7090 

EISP7091 

EISP7092 

EISP7093 

EISP7094 

EISP7095 

EISP7096 

EISP7097 

EISP7098 

EISP7099 
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EISP7105 

EISP7106 

EISP7107 

EISP7108 

EISP7109 

EISP7110 

EISP7111 
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260 CONTINUE 


GO TO 70 


END QZ STEP 


SE T ERROR — ALL EIGENVALUES HAVE NOT 
CONVERGED AFTER 30*N ITERATIONS • 

c ioo° ^ ERR .T. E ?. SAVE EPSB FOR USE BY QZVAL AND QZVEC . 
1001 IF (N .GT. 1) B (N , 1 ) = EPSB 

RETURN _ x 

C** THIS PROGRAM VALID ON FTN4 AND FTN5 ** 

END 

C ROUTINE NAME - PF262=QZVAL 

C FROM EISPACK 


C ROUTINE NAJ 

C FROM El SPA' 

C 

C 

C 

C LATEST REVISION 

C 

C 


- AUGUST 1,1984 

COMPUTER SCIENCES CORP. , HAMPTON, VA. 


EISP7112 
EISP7113 
EISP7114 
EISP7115 
EISP7116 
EISP7117 
EISP7118 
EISP7119 
EISP7120 
EISP7121 
EISP7122 
QZVAL 2 
QZVAL 3 
QZVAL 4 

QZVAL 5 

QZVAL 6 
QZVAL 7 
QZVAL 8 
QZVAL 9 
QZVAL 10 


C 

C PURPOSE 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C USAGE 
C 

C ARGUMENTS 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NM 

N 

A 


- THIS SUBROUTINE ACCEPTS A PAIR OF REAL 
MATRICES, ONE OF THEM IN QUASI -TRIANGULAR 
FORM AND THE OTHER IN UPPER TRIANGULAR FORM. 
IT REDUCES THE QUASI -TRIANGULAR MATRIX 
FURTHER, SO THAT ANY REMAINING 2-BY-2 BLOCKS 
CORRESPOND TO PAIRS OF COMPLEX EIGENVALUES, 
AND RETURNS QUANTITIES WHOSE RATIOS GIVE THE 
GENERALIZED EIGENVALUES. IT IS USUALLY 
PRECEDED BY QZHES(PF260) AND QZIT(PF261) AND 
MAY BE FOLLOWED BY QZVEC (PF2 63 ) . 


- CALL QZVAL (NM,N, A, B,ALFR,ALFI, BETA, MATZ,Z) 

- ON INPUT NM MUST BE SET TO THE ROW DIMENSION 
OF TWO-DIMENSIONAL ARRAY PARAMETERS AS 
DECLARED IN THE CALLING PROGRAM DIMENSION 
STATEMENT. 

- ON INPUT N IS THE ORDER OF THE MATRICES. 

- ON INPUT A CONTAINS A REAL UPPER QUASI- 
TRI ANGULAR MATRIX. 

MUST BE OF DIMENSION NM X N. 

ON OUTPUT A HAS BEEN REDUCED FURTHER TO A 
QUASI-TRIANGULAR MATRIX IN WHICH ALL NONZERO 
SUBDIAGONAL ELEMENTS CORRESPOND TO PAIRS OF 
COMPLEX EIGENVALUES. 

- ON INPUT B CONTAINS A REAL UPPER TRIANGULAR 
MATRIX. 

MUST BE OF DIMENSION NM X N . 

IN ADDITION, LOCATION B(N,1) CONTAINS THE 
TOLERANCE QUANTITY (EPSB) COMPUTED AND SAVED 
IN QZIT (PF261) . 


QZVAL 11 
QZVAL 12 
QZVAL 13 
QZVAL 14 
QZVAL 15 
QZVAL 16 
QZVAL 17 
QZVAL 18 
QZVAL 19 
QZVAL 20 
QZVAL 21 
QZVAL 22 
QZVAL 23 
QZVAL 24 
QZVAL 25 
QZVAL 26 
QZVAL 27 
QZVAL 28 
QZVAL 29 
QZVAL 30 
QZVAL 31 
QZVAL 32 
QZVAL 33 
QZVAL 34 
QZVAL 35 
QZVAL 36 
QZVAL 37 
QZVAL 38 
QZVAL 39 
QZVAL 40 
QZVAL 41 
QZVAL 42 
QZVAL 43 
QZVAL 44 
QZVAL 45 
QZVAL 46 
QZVAL 47 


ON OUTPUT B IS STILL IN UPPER TRIANGULAR QZVAL 48 
FORM, ALTHOUGH ITS ELEMENTS HAVE BEEN ALTERED . QZVAL 49 
B (N , 1 ) IS UNALTERED. QZVAL 50 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ALFR - ON OUTPUT ALFR CONTAINS THE REAL PART OF THE 
DIAGONAL ELEMENTS"" OF THE TRIANGULAR MATRIX 
THAT WOULD BE OBTAINED IF A WERE REDUCED 
COMPLETELY TO TRIANGULAR FORM BY UNITARY 
TRANSFORMATIONS. NON-ZERO VALUES OF ALFI 
OCCUR IN PAIRS, THE FIRST MEMBER POSITIVE AND 
THE SECOND NEGATIVE. 

MUST BE OF DIMENSION N. 

ALFI - ON OUTPUT ALFI CONTAINS THE IMAGINARY PART 


QZVAL 

QZVAL 

QZVAL 

QZVAL 

QZVAL 

QZVAL 

QZVAL 

QZVAL 

QZVAL 

QZVAL 

QZVAL 

QZVAL 


vnu 

OF THE DIAGONAL ELEMENTS OF OF THE TRIANGULAR QZVAL 


MATRIX THAT WOULD BE OBTAINED IF A WERE 
REDUCED COMPLETELY TO TRIANGULAR FORM BY 
UNITARY TRANSFORMATIONS. NON-ZERO VALUES 
OF ALFI OCCUR IN PAIRS, THE FIRST MEMBER 
POSITIVE AND THE SECOND NEGATIVE. 

MUST BE OF DIMENSION N. 


BETA 


QZVAL 
QZVAL 
QZVAL 
QZVAL 
QZVAL 
QZVAL 
QZVAL 

ON OUTPUT BETA CONTAINS THE DIAGONAL ELEMENTS QZVAL 
OF THE CORRESPONDING B, NORMALIZED TO BE REAL QZVAL 
AND NON-NEGATIVE. THE GENERALIZED EIGENVALUESQZVAL 


MATZ 


ARE THEN THE RATIOS ( (ALFR+I* ALFI) /BETA) . 
MUST BE OF DIMENSION N. 


- ON INPUT MATZ SHOULD BE SET TO .TRUE. IF 
THE RIGHT HAND TRANSFORMATIONS ARE TO BE 
ACCUMULATED FOR LATER USE IN COMPUTING 
EIGENVECTORS, AND TO .FALSE. OTHERWISE. 

Z - ON INPUT Z CONTAINS, IF MATZ HAS BEEN SET 

TO .TRUE., THE TRANSFORMATION MATRIX PRODUCED 
IN THE REDUCTIONS BY QZHES(PF260) AND QZIT 
(PF261) IF PERFORMED, OR ELSE THE IDENTITY 
MATRIX. IF MATZ HAS BEEN SET TO .FALSE. Z 
IS NOT REFERENCED. 

MUST BE OF DIMENSION NM X N. 

ON OUTPUT Z CONTAINS THE PRODUCT OF THE 
RIGHT HAND TRANSFORMATIONS (FOR ALL THREE 
STEPS) IF MATZ HAS BEEN SET TO .TRUE. 

REQUIRED ROUTINES - NONE 

REMARKS 1. THIS SUBROUTINE IS THE THIRD STEP OF THE QZ 
ALGORITHM FOR SOLVING GENERALIZED MATRIX 
EIGENVALUE PROBLEMS, SIAM J. NUMER. ANAL. 10, 
241-256(1973) BY MOLER AND STEWART. 

EXAMPLE : 

PROGRAM TQZVAL (OUTPUT, TAPE6=OUTPUT) 

LOGICAL°MATZ 5,5) ' B(5,5) ' ALFR(5) * ALFI ( 5 ) , BETA (5) , Z ( 5 , 5 ) 


51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 


N = 5 
NM = 5 
MATZ = 
EPS1 = 


. TRUE . 
0. 0E0 


QZVAL 
QZVAL 
QZVAL 
QZVAL 
QZVAL 
QZVAL 
QZVAL 
QZVAL 
QZVAL 82 
QZVAL 83 
QZVAL 84 
QZVAL 85 
QZVAL 86 
QZVAL 87 
QZVAL 88 
QZVAL 89 
QZVAL 90 
QZVAL 91 
QZVAL 92 
QZVAL 
QZVAL 
QZVAL 
QZVAL 
QZVAL 
QZVAL 
QZVAL 
QZVAL1 00 
QZVAL1 0 1 
QZVAL102 
QZVAL103 
QZVAL104 
QZVAL105 
QZVAL106 
QZVAL107 
QZVAL108 
QZVAL109 
QZVAL110 


93 

94 

95 

96 

97 

98 

99 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C99 

C100 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


DATA A /10.,2.,3.,2*1.,2.,12.,1.,2.,1.,3.,1.,11., 

DATA B /12. ,1. ,-l. ,2. ,2*1. i-1- '- 1 - i 1 - ' 

* 16. ,-l. ,1. ,2. ,-l. ,-l.,12. ,-i. ,3 *1- I -1 - r 11 • 

CALL QZHES (NM, N , A, B , MATZ , Z) 

CALL QZIT(NM,N,A,B,EPS1,MATZ,Z,IERR) 

CALL QZVAL (NM , N , A , B , ALFR , ALFI , BETA , MATZ , Z ) 

WRITE ( 6 \ 100) ALFR, ALFI , BETA, ( (Z(I,J) ,1=1,5) , J=l, 5 ) 
FORMAT ( 1H1 , 8H I ERR = ,14) 

FORMAT (1H , 8H ALFR = / 1H , 5(G8.2,2X)/ 

* 8H ALFI = /1H ,5(G8.2,2X)/ 

* 8H BETA = /1H , 5(G8.2,2X)/ 

* 5H Z = / 5 ( 1H , 5 (G8 . 2 , 2X) / ) ) 

STOP 

END 

OUTPUT : 


I ERR = 

0 




ALFR « 
15. 

7 . 2 

16. 

10. 

8.6 

ALFI = 
0 . 

0 . 

0 . 

0 . 

0 . 

BETA = 
9.9 

17. 

14. 

11. 

13. 

Z = 
.24 

- . 54E-01 

.21 

-.27 

-.91 


-.54 

.49 

-.60 

-.25 


.25 

.56 

.48 

-.63 


.65 

.49 

-.29 

.45 


-.46 

.45 

.44 

.57 


. 75E-01 
-.38 

- . 94E-01 


SUBROUTINE QZVAL ( NM , N , A , B , ALFR , ALFI , BETA , MATZ , Z) 

implicit real*8 (a-h,o-z) 

INTEGER I , J , N , EN , NA, NM, NN , ISW 

REAL* 8 A (NM, N) , B (NM, N) , ALFR (N) ALF! (N) BETA(N) ,Z ( NM, N) 

REAL* 8 C / D,E,R,S,T,AN,A1,A2,BH,CQ,CZ,DI,DR,EI,TI,TR,U1, 

X U2 , VI , V2 , All , All , A-12 , A2I , A21 , A22 , Bll , B12 , B22 , SQI , SQR , 

X SSI , SSR, SZI , SZR, A11I , A11R, A12I , A12R, A22I , A22R, EPSB 

LOGICAL MATZ 
EPSB = B (N , 1 ) 

T CJW = 1 

find EIGENVALUES OF QUASI -TRIANGULAR MATRICES. 

FOR EN=N STEP -1 UNTIL 1 DO 

DO 510 NN = 1, N 
EN = N + 1 - NN 
NA = EN - 1 

IF (ISW .EQ. 2) GO TO 505 

IF (EN .EQ. 1) GO TO 410 

IF (A(EN.NA) .NE. 0.0E0) GO TO 420 

1-BY-l BLOCK, ONE REAL ROOT 

410 " ’aLFR(EN) = A(EN,EN) 

IF (B (EN, EN) .LT. 0.0E0) ALFR(EN) - ALFR(EN) 

BETA (EN) = ABS (B (EN, EN) ) 


QZVAL111 
QZVAL112 
QZVAL113 
QZVAL114 
QZVAL115 
QZVAL116 
QZVAL117 
QZVAL118 
QZVAL119 
QZVAL120 
QZVAL12 1 
QZVAL122 
QZVAL123 
QZVAL124 
QZVAL125 
QZVAL126 
QZVAL127 
QZVAL128 
QZVAL129 
QZVAL130 
QZVAL131 
QZVAL132 
QZVAL133 
QZVAL134 
QZVAL135 
QZVAL136 
QZVAL137 
QZVAL138 
QZVAL139 
QZVAL140 
QZVAL14 1 
QZVAL142 
QZVAL143 
QZVAL144 
QZVAL145 
QZVAL146 
- QZVAL147 
EISP7123 

EISP7124 

EISP7125 

EISP7126 

EISP7127 

EISP7128 

EISP7129 

EISP7130 

EISP7131 

EISP7132 

EISP7133 

EISP7134 

EISP7135 

EISP713 6 

EISP7137 

EISP7138 

EISP7139 

EISP7140 

EISP7141 

EISP7142 

EISP7143 

EISP7144 
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o n 


0. OEO 


ALFI(EN) = 

GO TO 510 
2-BY-2 


BLOCK 


IF (ABS (B (NA, NA) ) . LE . EPSB) GO TO"” 455 
IF (ABS (B (EN , EN) ) . GT. EPSB) GO TO 430 
A1 = A ( EN , EN ) 

A2 = A (EN, NA) 

BN = 0.0E0 
GO TO 435 

AN = ABS (A (NA, NA) ) + ABS (A (NA , EN) ) + AB 
+ ABS (A (EN , EN) ) 

BN = ABS (B (NA, NA) ) + ABS (B (NA, EN) ) + AB 

All = A (NA, NA) / AN 

A12 = A (NA, EN) / AN 

A21 = A ( EN , NA) / AN 

A22 = A ( EN , EN ) / AN 

Bll = B (NA, NA) / BN 

B12 = B (NA, EN) / BN 

B22 = B ( EN , EN ) / BN 

E = All / Bll 

El = A22 / B22 

S = A21 / (Bll * B22) 

T = (A22 - E * B22 ) / B22 

IF (ABS (E) .LE. ABS ( El ) ) GO TO 431 

E = El 

T = (All - E * Bll) / Bll 
C = 0.5E0 * (T - S * B12 ) 

D = C*C + S* ( A12 - E * B12 ) 

IF (D .LT. 0.0E0) GO TO 480 

■ TWO REAL ROOTS . 

ZERO BOTH A ( EN , NA) AND B ( EN , NA) 
E = E + (C + SIGN ( SQRT ( D) , C) ) 


GO TO"” 455 
GO TO 430 


+ ABS (A (NA , EN) ) + ABS (A (EN, NA) ) 
+ ABS (B (NA, EN) ) + ABS (B ( EN , EN) ) 


* B22) 

B22 ) / B22 
. ABS (El) ) GO TO 431 


D = C 
IF (D 


All 

= All - 

E 

* Bll 

A12 

= A12 - 

E 

* B12 

A22 

= A22 - 

E 

* B22 

IF 

(ABS (All) 

-f 

ABS (A12 ) 


ABS (A2 1) 
= A12 
= All 


ABS (A12 ) .LT. 

ABS (A22 ) ) GO TO 432 


GO TO 


= A22 
= A21 


CHOOSE AND APPLY REAL 

S = ABS(Al) + ABS ( A2 ) 

U1 = A1 / S 
U2 = A2 / S 

R = SIGN(SQRT(U1*U1+U2*U2) ,U1) 
VI = -(Ul + R) / R 
V2 = -U2 / R 
U2 = V2 / VI 


DO 440 I = 1, EN 

T = A ( I , EN ) + U2 * A ( I , NA) 
A ( I , EN) = A ( I , EN) + T * VI 

A ( I , NA) = A(I,NA) 4- T * V2 

T = B (I , EN) + U2 * B ( I , NA) 
B ( I , EN) = B ( I , EN ) + T * VI 

B ( I , NA) = B ( I , NA) + T * V2 

CONTINUE 


EISP7145 

EISP7146 

EISP7147 

EISP7148 

EISP7149 

EISP7150 

EISP7151 

EISP7152 

EISP7153 

EISP7154 

EISP7155 

EISP7156 

EISP7157 

EISP7158 

EISP7159 

EISP7160 

EISP7161 

EISP7162 

EISP7163 

EISP7164 

EISP7165 

EISP7166 

EISP7167 

EISP7168 

EISP7169 

EISP7170 

EISP7171 

EISP7172 

EISP7173 

EISP7174 

EISP7175 

EISP7176 

EISP7177 

EISP7178 

EISP7179 

EISP7180 

EISP7181 

EISP7182 

EISP7183 

EISP7184 

EISP7185 

EISP7186 

EISP7187 

EISP7188 

EISP7189 

EISP7190 

EISP7191 

EISP7192 

EISP7193 

EISP7194 

EISP7195 

EISP7196 

EISP7197 

EISP7198 

EISP7199 

EISP7200 

EISP72 01 

EISP7202 

EISP7203 

EISP7204 
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IF (.NOT. MATZ) GO TO 450 
C 

DO 445 I = 1, N 

T = Z ( I , EN ) + U2 * Z ( I , NA) 

Z ( I , EN ) = Z ( I , EN) + T * VI 
Z ( I , NA ) = Z ( I , NA) + T * V2 
445 CONTINUE 

C 

450 IF (BN .EQ. 0.0E0) GO TO 475 

IF (AN .LT. ABS(E) * BN) GO TO 455 

A1 = B (NA, NA) 

A2 = B (EN, NA) 

GO TO 460 

455 A1 = A (NA, NA) 

A2 = A (EN, NA) 

c CHOOSE AND APPLY REAL Q 

460 " S = ABS(Al) + ABS ( A2 ) 

IF (S .EQ. O.OEO) GO TO 475 
U1 = A1 / S 

U2 = A2 / S 

R = SIGN (SQRT (U1*U1+U2*U2 ) ,U1) 

VI = -(Ul + R) / R 
V2 = -U2 / R 

U2 = V2 / VI 

C 

DO 470 J = NA, N 

T = A(NA, J) + U2 * A ( EN , J ) 

A(NA, J) = A (NA, J) + T * VI 

A ( EN , J ) = A ( EN , J ) + T * V2 

T = B (NA , J) + U2 * B ( EN , J ) 

B (NA, J) = B (NA, J) + T * VI 

B ( EN , J) = B ( EN , J ) + T * V2 

470 CONTINUE 

C 

475 A(EN, NA) = O.OEO 

B(EN,NA) = O.OEO 

ALFR(NA) = A(NA,NA) 

ALFR(EN) = A (EN , EN) 

IF (B (NA, NA) .LT. O.OEO) ALFR(NA) = -ALFR(NA) 

IF ( B ( EN , EN ) .LT. O.OEO) ALFR(EN) = -ALFR(EN) 

BETA (NA) = ABS (B (NA, NA) ) 

BETA(EN) = ABS (B (EN, EN) ) 

ALFI (EN) = O.OEO 
ALFI (NA) = O.OEO 
GO TO 505 

C TWO COMPLEX ROOTS 

480 E = E + C 

El = SQRT(-D) 

A11R = All - E * Bll 
A11I = El * Bll 
A12R = A12 - E * B12 
A12I = El * B12 
A22R = A22 - E * B22 
A22I = El * B22 

IF (ABS(AllR) + ABS(AllI) + ABS(A12R) + ABS(A12I) .LT. 
X ABS ( A2 1) + ABS (A22R) + ABS(A22I)) GO TO 482 

A1 = A12R 
All = A12I 
A2 = -A11R 
A2I = -A11I 


EISP7205 

EISP7206 

EISP7207 

EISP7208 

EISP7209 

EISP7210 

EISP7211 

EISP7212 

EISP72 13 

EISP7214 

EISP72 15 

EISP7216 

EISP7217 

EISP7218 

EISP7219 

EISP7220 

EISP7221 

EISP7222 

EISP7223 

EISP7224 

EISP7225 

EISP7226 

EISP7227 

EISP7228 

EISP7229 

EISP7230 

EISP7231 

EISP7232 

EISP7233 

EISP7234 

EISP7235 

EISP7236 

EISP7237 

EISP7238 

EISP7239 

EISP7240 

EISP7241 

EISP7242 

EISP7243 

EISP7244 

EISP724 5 

EISP7246 

EISP7247 

EISP7248 

EISP7249 

EISP7250 

EISP7251 

EISP7252 

EISP7253 

EISP7254 

EISP7255 

EISP7256 

EISP7257 

EISP7258 

EISP7259 

EISP7260 

EISP7261 

EISP7262 

EISP7263 

EISP7264 
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o n 


482 


GO TO 485 
A1 = A22R 
All = A22I 
A2 = -A21 
A2I = 0.0E0 

C CHOOSE COMPLEX Z 

485 CZ = SQRT (A1*A1+A1I*A1I) 

IF (CZ .EQ. 0.0E0) GO TO 487 
SZR = ( A1 * A2 + All * A2I ) / CZ 

SZI = ( A1 * A2I - All * A2 ) / CZ 

R = SQRT (CZ*CZ+SZR*SZR+SZI*SZI) 
CZ = CZ / R 
SZR = SZR / R 

SZI = SZI / R 


GO TO 490 
487 SZR = 1.0E0 
SZI = O.OEO 

490 IF (AN .LT. (ABS(E) + El) * BN) GO TO 492 
A1 = CZ * Bll + SZR * B12 
All = SZI * B12 
A2 = SZR * B22 
A2I = SZI * B22 


GO TO 495 

492 A1 = CZ * All + SZR * A12 
All = SZI * A12 
A2 = CZ * A21 + SZR * A22 
A2I = SZI * A22 

C CHOOSE COMPLEX Q 

495 CQ = SQRT (A1*A1+A1I*A1I ) 

IF (CQ .EQ. O.OEO) GO TO 497 
SQR = (A1 * A2 + All * A2I ) / CQ 

SQI = ( A1 * A2I - All * A2 ) / CQ 

R = SQRT ( CQ*CQ+SQR*SQR+SQI *SQI ) 
CQ = CQ / R 
SQR = SQR / R 

SQI = SQI / R 


GO TO 500 

497 SQR = l.OEO 

SQI = O.OEO 

COMPUTE DIAGONAL ELEMENTS THAT WOULD RESULT 

IF TRANSFORMATIONS WERE APPLIED 

500 SSR = SQR * SZR + SQI * SZI 

SSI = SQR * SZI - SQI * SZR 

1 = 1 

TR = CQ * CZ * All + CQ * SZR * A12 + SQR * CZ * A21 

X + SSR * A22 

TI = CQ * SZI * A12 - SQI * CZ * A21 + SSI * A22 

DR = CQ * CZ * Bll + CQ * SZR * B12 + SSR * B22 

DI = CQ * SZI * B12 + SSI * B22 

GO TO 503 

502 1=2 

TR = SSR * All - SQR * CZ * A12 - CQ * SZR * A21 
X + CQ * CZ * A22 

TI = -SSI * All - SQI * CZ * A12 + CQ * SZI * A21 
DR = SSR * Bll - SQR * CZ * B12 + CQ * CZ * B22 
DI = -SSI * Bll - SQI * CZ * B12 

503 T = TI * DR - TR * DI 
J = NA 

IF (T .LT. O.OEO) J = EN 
R = SQRT (DR*DR+DI*DI ) 


EISP7265 

EISP7266 

EISP72 67 

EISP7268 

EISP7269 

EISP7270 

EISP7271 

EISP7272 

EISP7273 

EISP7274 

EISP7275 

EISP7276 

EISP7277 

EISP7278 

EISP7279 

EISP7280 

EISP7281 

EISP7282 

EISP7283 

EISP7284 

EISP7285 

EISP7286 

EISP7287 

EISP7288 

EISP7289 

EISP7290 

EISP7291 

EISP7292 

EISP7293 

EISP7294 

EISP7295 

EISP7296 

EISP7297 

EISP7298 

EISP7299 

EISP7300 

EISP73 01 

EISP7302 

EISP7303 

EISP7304 

EISP7305 

EISP7306 

EISP7307 

EISP7308 

EISP7309 

EISP7310 

EISP73 11 

EISP73 12 

EISP7313 

EISP7314 

EISP7315 

EISP7316 

EISP7317 

EISP7318 

EISP7319 

EISP7320 

EISP7321 

EISP7322 

EISP7323 

EISP7324 
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LATEST REVISION 


PURPOSE 


BETA ( J) = BN * R 

ALFR(J) = AN * (TR * DR + TI * DI) / R 
ALFI(J) = AN * T / R 
IF (I .EQ. 1) GO TO 502 

505 ISW = 3 - ISW 
510 CONTINUE 

B (N, 1) = EPSB 
C 

RETURN 

C** THIS PROGRAM VALID ON FTN4 AND FTN5 ** 

END 

C ROUTINE NAME - PF263=QZVEC 

C FROM EISPACK 

C 



C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


USAGE 

ARGUMENTS NM 


N 

A 


EISP7325 
EISP7326 
EISP7327 
EISP7328 
EISP7329 
EISP7330 
EISP73 3 1 
EISP7332 
EISP7333 
EISP7334 
EISP73 35 
QZVEC 2 
QZVEC 3 
QZVEC 4 
QZVEC 5 


- AUGUST 1,1984 

COMPUTER SCIENCES CORP. , HAMPTON, VA. 

- THIS SUBROUTINE ACCEPTS A PAIR OF REAL 
MATRICES, ONE OF THEM IN QUASI -TRIANGULAR 
FORM (IN WHICH EACH 2-BY-2 BLOCK CORRESPONDS 
TO A PAIR OF COMPLEX EIGENVALUES) AND THE 
OTHER IN UPPER TRIANGULAR FORM. IT COMPUTES 
THE EIGENVECTORS OF THE TRIANGULAR PROBLEM 
AND TRANSFORMS THE RESULTS BACK TO THE 
ORIGINAL COORDINATE SYSTEM. IT IS USUALLY 
PRECEDED BY QZHES (PF260) , QZIT(PF261), AND 
QZVAL (PF2 62 ) . 

- CALL QZVEC (NM , N , A, B , ALFR, ALFI , BETA, Z) 

- ON INPUT NM MUST BE SET TO THE ROW DIMENSION 
OF TWO-DIMENSIONAL ARRAY PARAMETERS AS 
DECLARED IN THE CALLING PROGRAM DIMENSION 
STATEMENT . 

- ON INPUT N IS THE ORDER OF THE MATRICES. 

- ON INPUT A CONTAINS A REAL UPPER QUASI- 
TRI ANGULAR MATRIX. 

MUST BE OF DIMENSION NM X N . 


6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 


QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 17 
QZVEC 18 
QZVEC 19 
QZVEC 20 
QZVEC 21 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 28 
QZVEC 29 
QZVEC 30 
QZVEC 31 
QZVEC 32 
QZVEC 33 
QZVEC 34 
QZVEC 35 
QZVEC 36 


22 

23 

24 

25 

26 
27 


ON OUTPUT A IS UNALTERED. ITS SUBDIAGONAL 
ELEMENTS PROVIDE INFORMATION ABOUT THE STORAGEQZVEC 37 

OF THE COMPLEX EIGENVECTORS. QZVEC 38 

QZVEC 39 

ON INPUT B CONTAINS A REAL UPPER TRIANGULAR QZVEC 40 
MATRIX. IN ADDITION, LOCATION B(N,1) CONTAINSQZVEC 41 
THE TOLERANCE QUANTITY (EPSB) COMPUTED AND 
SAVED IN QZIT (PF261 ) . 

MUST BE OF DIMENSION NM X N. 


ALFR 


ON OUTPUT B HAS BEEN DESTROYED. 

ON INPUT ALFR IS A VECTOR SUCH THAT THE 
RATIOS ( (ALFR+I*ALFI) /BETA) ARE THE 
GENERALIZED EIGENVALUES. THEY ARE USUALLY 


QZVEC 42 
QZVEC 43 
QZVEC 44 
QZVEC 45 
QZVEC 46 
QZVEC 47 
QZVEC 48 
QZVEC 49 
QZVEC 50 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ALFI 


BETA 


REQUIRED ROUTINES 
REMARKS 1 


OBTAINED FROM QZVAL(PF262) . QZVEC 51 

MUST BE OF DIMENSION N. QZVEC 52 

QZVEC 53 

ON INPUT ALFI IS^ VECTOR SUCH THAT THE RATIOSQZVEC 54 
( (ALFR+I*ALFI) /BETA) ARE THE GENERALIZED QZVEC 55 

EIGENVALUES. THEY ARE USUALLY OBTAINED FROM QZVEC 56 
QZVAL(PF262) . QZVEC 57 

MUST BE OF DIMENSION N. QZVEC 58 

QZVEC 59 

- ON INPUT BETA IS A VECTOR_SUCH_THAT THE RATIOSQZVEC 60 

QZVEC 61 
QZVEC 62 
QZVEC 63 
QZVEC 64 
QZVEC 6 5 

ON INPUT Z CONTAINS THE TRANSFORMATION MATRIX QZVEC 66 
PRODUCED IN THE REDUCTIONS BY QZHES (PF260) QZVEC 67 
QZIT(PF261), AND QZVAL(PF262) , IF PERFORMED. QZVEC 68 
IF THE EIGENVECTORS OF THE TRIANGULAR PROBLEM QZVEC 69 


A V1\ u UUU ± ILn X InL i 

( (ALFR+I*ALFI) /BETA) ARE THE GENERALIZED 
EIGENVALUES. THEY ARE USUALLY OBTAINED FROM 
QZVAL ( PF2 62 ) . 

MUST BE OF DIMENSION N. 


ARE DESIRED, Z MUST CONTAIN THE IDENTITY 
MATRIX. 

MUST BE OF DIMENSION NM X N. 


QZVEC 70 
QZVEC 71 
QZVEC 72 
QZVEC 73 
QZVEC 74 
QZVEC 75 


ON OUTPUT Z CONTAINS THE REAL AND IMAGINARY 
PARTS OF THE EIGENVECTORS. IF ALFI (I) .EQ. 

0.0, THE I-TH EIGENVALUE IS REAL AND THE I-TH QZVEC 76 
COLUMN OF Z CONTAINS ITS EIGENVECTOR. IF QZVEC 77 

ALFI (I) .NE. 0.0, THE I-TH EIGENVALUE IS QZVEC 78 

to M ™« EX * IF ALFI ( 1 ) - GT * 0.0, THE EIGENVALUE QZVEC 79 
J® FIRST 0F A COMPLEX PAIR AND THE I-TH QZVEC 80 

AND (I+1)-TH COLUMNS OF Z CONTAIN ITS EIGEN- QZVEC 81 
VECTOR. IF ALFI (I) .LT. 0.0, THE EIGEN- QZVEC 82 

VALUE IS THE SECOND OF A COMPLEX PAIR AND THEQZVEC 83 

QZVEC 84 


( I- 1)~TH AND I-TH COLUMNS OF Z CONTAIN THE 
CONJUGATE OF ITS EIGENVECTOR. EACH EIGEN- 
VECTOR IS NORMALIZED SO THAT THE MODULUS 
OF ITS LARGEST COMPONENT IS 1.0 . 


- NONE 


THIS SUBROUTINE IS THE OPTIONAL FOURTH STEP 
OF THE QZ ALGORITHM FOR SOLVING GENERALIZED 
MATRIX EIGENVALUE PROBLEMS, SIAM J. NUMER. 
ANAL. 10, 241-256(1973) BY MOLER AND STEWART. 

EXAMPLE : 

PROGRAM TQZVEC ( OUTPUT , TAPE6=OUTPUT) 

LOGICAL°MATZ 5,5) ,B(5,5) ,ALFR(5) ' ALFI (5) , BETA (5) , Z (5 , 5 ) 


QZVEC 85 
QZVEC 86 
QZVEC 87 
QZVEC 88 
QZVEC 89 
QZVEC 90 
QZVEC 91 
QZVEC 92 
QZVEC 93 
QZVEC 94 
QZVEC 95 
QZVEC 96 
QZVEC 97 
QZVEC 98 
QZVEC 99 
QZVEC100 


c 

N = 5 


QZVEC101 

c 

NM = 5 


QZVEC102 

c 

MATZ = 

. TRUE . 

QZVEC103 

c 

EPS1 = 

0. 0E0 

QZVEC104 

c 



QZVEC105 

c 

c 

c 

DATA A 

* 

/ 10 . , 2 . , 3 . ,2*1. ,2. ,12. ,1. ,2. ,1. , 3 . ,1. ,11. , 
1. ,-l. , 1. , 2. , 1. ,9. ,3*1. ,-l. , l. , 15. / 

QZVEC106 

QZVEC107 

QZVEC108 

c 

DATA B 

/12. ,1. , -1 . ,2. ,2*1. ,14. ,1. , -1 . ,1. ,-l. ,1. t 

QZVEC109 

QZVEC110 
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C * 16. ,-l. ,1. , 2. ,-l. ,-l- / 12 - I -1 * f 3 * 1 - *" 1 ‘ ' 11 ‘ 1 

C CALL QZHES ( NM , N , A , B , MATZ , Z ) 

C CALL QZIT(NM,H,A,B,EPSl,MATZ,Z,IERRr“ 

C CALL QZVAL(NM, N , A, B , ALFR, ALFI , BETA, MATZ , Z) 

C CA LL QZVEC(NM, N, A, B, ALFR, ALFI, BETA, Z) 

C WRITE (6,99) I ERR 

C WRITE ( 6 , 100 ) ( (Z(I,J) ,1=1,5) ,J-1,5) 

C99 FORMAT ( 1H1 , 7HIERR = ,14) 

C100 FORMAT ( 5H Z = /5(1H , 5 (G8 . 2 , 2X) / ) ) 

C STOP 

C END 

C 

C OUTPUT : 

C 

C I ERR = 0 

C %6 -- 59E-01 .23 -30 -1.0 

l SI i”. ^ -• *H E-01 

\ -:S -:S i!i 

SUBROUTINE QZVEC (NM , N , A, B, ALFR, ALFI , BETA, Z) 

C 

INTEGER^i r j a K,M,N, EN* II , JJ, NA, NM, NN , ISW, ENM2 

REAL* 8 A(Nm',N) ,B(NM,N) , ALFR(N), ALFI (N) , BETA ( N^ , Z(NM^ 
REAL* 8 D,Q,R,S,T,W,X,Y,DI,DR,RA,RR,SA,TI,TR,T1,T2,W1,X1, 

X Z Z , Z 1 , ALFM , ALMI , ALMR , BETM , EPSB 

EPSB = B (N , 1) 

XSW = 1 

c for EN=N STEP -1 UNTIL 1 DO — 

DO 800 NN = 1, N 
EN = N + 1 - NN 
NA = EN - 1 

IF (ISW .EQ. 2) GO TO 795 
IF (ALFI (EN) .NE. 0.0E0) GO TO 710 

c REAL VECTOR 

M = EN 

B ( EN , EN ) = 1.0E0 

IF (NA .EQ. 0) GO TO 800 

ALFM = ALFR (M) 

BETM = BETA (M) 

c FOR I=EN-1 STEP -1 UNTIL 1 DO — 

DO 700 II = 1, NA 
I = EN - II 

W = BETM * A (I, I) “ ALFM * B(I,I) 

R = 0.0E0 

610 r°= 6 R 0 + J (BETM * N A(I,J) ~ ALFM * B(I,J)) * B(J,EN) 

C IF (I .EQ. 1 -OR. ISW .EQ. 2) GO TO 630 

IF (BETM * A(I,I~1) .EQ. 0.0E0) GO TO 630 
ZZ = W 
S = R 
GO TO 690 
630 M = I 


QZVEC1 1 1 

QZVEC112 

QZVEC113 

QZVEC114 

QZVEC115 

QZVEC116 

QZVEC117 

QZVEC118 

QZVEC119 

QZVEC120 

QZVEC121 

QZVEC122 

QZVEC123 

QZVEC124 

QZVEC125 

QZVEC126 

QZVEC127 

QZVEC128 

QZVEC129 

QZVEC130 

QZVEC131 

QZVEC132 

QZVEC133 

QZVEC134 

EISP7336 

EISP7337 

EISP7338 

EISP7339 

EISP7340 

EISP734 1 

EISP7342 

EISP7343 

EISP7344 

EISP7345 

EISP7346 

EISP7347 

EISP7348 

EISP7349 

EISP7350 

EISP7351 

EISP7352 

EISP7353 

EISP7354 

EISP7355 

EISP7356 

EISP7357 

EISP7358 

EISP7359 

EISP7360 

EISP7361 

EISP7362 

EISP7363 

EISP7364 

EISP7365 

EISP7366 

E1SP7367 

EISP7368 

EISP7369 

EISP7370 
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o o 


IF (ISW .EQ. 2) GO TO 640 

C REAL 1-BY-l BLOCK 

T = W 

IF (W .EQ. 0.0E0) T = EPSB 
B ( I , EN ) = -R / T 
GO TO 700 

C REAL 2-BY-2 BLOCK 

640 X = BETM * A(I,I+1) - ALFM * B(I,I+1) 

Y = BETM * A (1+1,1) 


Q = W * ZZ - X * Y 
T = (X * S - ZZ * R) / Q 
B ( I , EN) = T 

IF (ABS(X) .LE. ABS(ZZ)) GO TO 650 
B ( 1+1 , EN) = (-R - W * T) / X 
GO TO 690 


650 B ( 1+1 , EN) = (-S - Y * T) / ZZ 

690 ISW = 3 - ISW 

700 CONTINUE 

C END REAL VECTOR 

GO TO 800 

C COMPLEX VECTOR 

710 M = NA 


ALMR = ALFR(M) 

ALMI = ALFI(M) 

BETM = BETA (M) 

LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT 

EIGENVECTOR MATRIX IS TRIANGULAR 

Y = BETM * A (EN, NA) 

B (NA, NA) = -ALMI * B ( EN , EN ) / Y 

B (NA, EN) = (ALMR * B ( EN , EN ) - BETM * A (EN, EN) ) / Y 
B(EN,NA) = 0.0E0 
B ( EN , EN ) = 1.0E0 
ENM2 = NA - 1 

IF (ENM2 .EQ. 0) GO TO 795 

FOR I=EN-2 STEP -1 UNTIL 1 DO — 

DO 790 II = 1, ENM2 
I = NA - II 

W = BETM * A (I, I) - ALMR * B(I,I) 

W1 = -ALMI * B(I,I) 

RA = 0.0E0 
SA = O.OEO 

DO 760 J = M, EN 

X = BETM * A(I,J) - ALMR * B(I,J) 

XI = -ALMI * B ( I , J) 

RA = RA + X * B ( J, NA) - XI * B(J,EN) 

SA = SA + X * B ( J , EN) + XI * B(J,NA) 

760 CONTINUE 

C 

IF (I .EQ. 1 .OR. ISW .EQ. 2) GO TO 770 
IF (BETM * A ( I , 1-1 ) .EQ. O.OEO) GO TO 770 
ZZ = W 
Z1 = W1 
R = RA 
S = SA 
ISW = 2 
GO TO 790 


770 M = I 

IF (ISW .EQ. 2) GO TO 780 
c COMPLEX 1-BY-l BLOCK 


EISP7371 

EISP7372 

EISP7373 

EISP7374 

EISP7375 

EISP7376 

EISP7377 

EISP7378 

EISP7379 

EISP7380 

EISP7381 

EISP7382 

EISP7383 

EISP7384 

EISP7385 

EISP7386 

EISP7387 

EISP7388 

EISP7389 

EISP7390 

EISP7391 

EISP7392 

EISP7393 

EISP7394 

EISP7395 

EISP7396 

EISP7397 

EISP7398 

EISP7399 

EISP7400 

EISP7401 

EISP7402 

EISP7403 

EISP7404 

EISP7405 

EISP7406 

EISP7407 

EISP7408 

EISP7409 

EISP74 10 

EISP7411 

EISP7412 

EISP7413 

EISP74 14 

EISP7415 

EISP7416 

EISP7417 

EISP74 18 

EISP7419 

EISP7420 

EISP7421 

EISP7422 

EISP7423 

EISP7424 

EISP7425 

EISP7426 

EISP7427 

EISP7428 

EISP7429 

EISP7430 
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non o o o non 


773 


TR = -RA 
TI = -SA 
DR = W 
DI = W1 


C COMPLEX DIVIDE (T1,T2) = (TR,TI) / (DR, DI) 

775 IF (ABS (DI) .GT. ABS (DR) ) GO TO 777 


RR = DI / DR 
D = DR + DI * RR 
Tl = (TR + TI * RR) / D 

T2 = (TI - TR * RR) / D 

GO TO (787,782) , ISW 
CALL GOTOER 
777 RR = DR / DI 

D = DR * RR + DI 
Tl = (TR * RR + TI) / D 

T2 = (TI * RR - TR) / D 

GO TO (787,782) , ISW 
CALL GOTOER 


C COMPLEX 2-BY-2 BLOCK 

780 X = BETM * A(I,I+1) - ALMR * B(I,I+1) 

XI = -ALMI * B (I , 1+1) 

Y = BETM * A ( 1+1 , I ) 


TR = Y * RA -W * R+Wl * S 

TI = Y * SA -W*S-W1*R 

DR = W * ZZ - W1 * Z1 - X * Y 

DI = W * Z1 + W1 * ZZ - XI * Y 

IF (DR .EQ. 0.0E0 .AND. DI .EQ. 0.0E0) DR = EPSB 
GO TO 775 

782 B ( 1+1 , NA) = Tl 

B ( 1+1 , EN) = T2 
ISW = 1 

IF (ABS ( Y) .GT. ABS ( W) + ABS(Wl)) GO TO 785 
TR = -RA - X * B ( 1+1 , NA) + XI * B(I+1,EN) 

TI = -SA - X * B ( I + 1 , EN ) - XI * B ( 1+1 , NA) 

GO TO 773 


785 Tl = ( -R - ZZ * B(I+1,NA) + Z1 * B(I+1,EN)) / Y 

T2 = ( -S - ZZ * B ( 1+1 , EN) - Z1 * B(I+1,NA)) / Y 
787 B ( I , NA) = Tl 

B ( I , EN) = T2 
790 CONTINUE 

C END COMPLEX VECTOR 

795 ISW = 3 - ISW 


800 CONTINUE 

END BACK SUBSTITUTION. 

TRANSFORM TO ORIGINAL COORDINATE SYSTEM. 

FOR J=N STEP -1 UNTIL 1 DO — 

DO 880 JJ = 1, N 
J = N + 1 - JJ 

DO 880 I = 1, N 
ZZ = 0.0E0 

DO 860 K = 1, J 

860 ZZ = ZZ + Z ( I , K) * B (K , J) 

Z(I,J) = ZZ 
880 CONTINUE 

NORMALIZE SO THAT MODULUS OF LARGEST 

COMPONENT OF EACH VECTOR IS 1. 

(ISW IS 1 INITIALLY FROM BEFORE) 


EISP7431 

EISP7432 

EISP7433 

EISP7434 

EISP7435 

EISP7436 

EISP7437 

EISP7438 

EISP7439 

EISP7440 

EISP7441 

EISP7442 

EISP7443 

EISP7444 

EISP7445 

EISP7446 

EISP7447 

EISP7448 

EISP7449 

EISP7450 

EISP7451 

EISP7452 

EISP7453 

EISP7454 

EISP7455 

EISP7456 

EISP7457 

EISP7458 

EISP7459 

EISP7460 

EISP7461 

EISP7462 

EISP7463 

EISP7464 

EISP7465 

EISP7466 

EISP7467 

EISP7468 

EISP7469 

EISP7470 

EISP7471 

EISP7472 

EISP7473 

EISP7474 

EISP7475 

EISP7476 

EISP7477 

EISP7478 

EISP7479 

EISP7480 

EISP7481 

EISP7482 

EISP7483 

EISP7484 

EISP7485 

EISP7486 

EISP7487 

EISP7488 

EISP7489 

EISP7490 
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ouo6 


c 

c 


890 


900 


920 


930 


C 

C 


00 ”0 J - 1, H 
D = 0 . 0 E 0 

5 {Stiff- *> » » «. 

' * NE * 0 • OEO) GO TO 945 

DO 890 I = lf N 

CONTINUE ^ Z ( I ' J )) ■ GT . D) D == ABS (2 (I , j) j 
DO 900 I=i N 

- z(i;j? / D 

GO TO 950 
DO 930 I = lf N 

IF~(f S ^ a 6foEoj R “ 

CONTINUE * GT - D) ° = R + ^d,J)/R)i* 2) 

DO 940 I = 1 , N 

Z(t't7 1) = Z ( I >J-1) / D 

CONTINUE ‘ 2(I ' J) / D 

I?" * 3 - ISW 


940 

ion = 

950 CONTINUE 

RETURN 

END 

ROUTINE NAME - pr, c , _ 

FROM EISPACK PF266=Rgg 


EISP7491 

EISP7492 

EISP7493 

EISP7494 

EISP7495 

EISP7496 

EISP7497 

EISP7498 

EISP74 99 

EISP7500 

EISP7501 

EISP7502 

EISP7503 

EISP 7504 

EISP7505 

EISP7506 

EISP7507 

EISP7508 

EISP7509 

EISP7510 

EISP7511 

EISP7512 

EISP7513 

EISP7514 

EISP75I5 

EISP7516 

EISP7517 

EISP7518 

EISP7519 

EISP7520 

EISP7521 

"Ti 


c 

c 

c 

c 

c 

c 

c 


latest revision 

PURPOSE 


USAGE 

ARGUMENTS Njj 


N 


- AUGUST 1,1984 

COMPUTER SCIENCES CORE 

CUKP - ' HAMPTON, VA. 

THIS SUBROUTINP put. 

sequence op sutoSS^Sp™? recommended 
“I, INE PACKAGE (EiIpS ™ e eic ensistem 
the Es and eigenvectors T?„ P ™ D ™ e 

- ‘JSS&SS^ G “™«?E R D S ttS aT 

---Rcc ( n„,n,a.b.auer, ALF i. BETAiMatZi2i 
OF ™e U tw™dwensional T a TO ™ e eow dimension 

- « XNPUT n IS TME order OP the matrpces a 
' Wo? SS^Cn”^ matrix. 


RGG 
RGG 
RGG 
■ RGG 
RGG 
RGG 
RGG 
RGG 


2 

3 

4 

5 

6 

7 

8 

9 


RGG 10 
RGG 11 
RGG 12 
RGG 13 
RGG 1 4 
RGG 15 
RGG 16 
RGG 17 
RGG 18 
RGG 1 9 
RGG 20 
RGG 21 
RGG 22 
RGG 23 
RGG 24 
RGG 25 
RGG 2 6 
RGG 27 
RGG 28 
RGG 29 
RGG 30 
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o o 



B 


_ _ rvNERAL MATRIX* 

' SSot’m'o? SKsIOH^HM X »• 

-ritn ffU F 


ALFR 


' JSsi BE OF DIME*— 

' oN 0 UT tSrs^f R the N ^ ^ige^valSes. 
numerators d^ nsion n. 

MUST BE OF uj - -nfr n 


RGG 

RGG 

RGG 

RGG 


RGG 

RGG 


31 

32 

33 

34 

35 

36 

37 


ALFI 


BETA 


matz 


ox ° OTFOT gs srSFaSwiS**- 

the numerators OT sj TH n _ 

MUST s the denominators 

OH OUTPUT BETA COKTAWS^T ^ GIVE N 

™ E THE G rTIOS (AEFRPI'AEFII /““^es API 

Se^bS“«™ ^oinSfpwT first. 

^TIteoer variaeee set 
on INPUT MATZ IS ^ geNVAL UES ARE 
TO ZERO IF ONLY SE T TO 

«P22:„S , 3S£r for both ETC 


RGG 

38 

RGG 

39 

RGG 

40 

RGG 

41 

RGG 

42 

RGG 

43 

RGG 

44 

RGG 

45 

’ RGG 

46 

RGG 

47 

RGG 

48 

RGG 

49 

RGG 

50 

alrgg 

51 

c. o 

RGG 

52 

RGG 

53 


ON lNrux - v KlGENVAL.ur,^ 

eicehvalues aho - 
ANY NON -ZERO IHTW»«» r G G 

EIGENVECTORS. XMAGINARY rug 

rpur REAL AND iron r G G 


BD ROUTINES 


... RGG 

EIGENVECTORS. XMAGINARY P^G 

ON OUTPUT HOG 

is* "JvHsss-s! >* *“• s 

J-TH COLUMN OF Z ^ j_ th ve xmAGINARY RGG 

eigenvector. £ gmplex with posit RG g 

EIGENVALUE I (J+1 )-TH RGG 

PART, THE th e real and rgg 

COLU^S OF ITS EIGENVECTOR. THE R GG 

IMAGINARY PAR VECTOR IS TH p TGEN VALUE . RGG 

CONJUGATE OF THIS coNJU GATE EIGENVA RCC 

EIGENVECTOR FOR T nM X N. RGG 

otst be of uimbh ^ er ompuT ble POO 

_ ON OUTPUT IBR I » COMPLETION COD R G G 

SET EQUAL TO A» DOCUMENTATION FOR Z£R0> rGG 

?f 2 61) BE ?HE FORMAL COMPLETION * GG 

fil 0ZIT PF262=QZVAL,PF263 =QZVECRG g G g 
_ p F 260=QZHES,PF261=QZIT,P rgg 

HC318=EPSL0N rg G 

RGG 

references eige hsyste« ROUTIHES. boo 

EROH THE El SPACE PACKA ^ CBLL S ROUTINES^ 

-sissr z — «■ ■ « 

OZVEC (PF263 ) . rPNERAL MATRI CESRGG 


55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 


I ERR 


69 

70 

71 

72 

73 
7 A 
75 
71 
T 
1 


1 

8 

8 

* 


172 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

C99 

C100 


RGG 

Q2IT (PF2 61) ACCEPTS A PAIR OF REAL MATRICES, ONE OFRGG 

THEM IN UPPER HESSENBERG FORM AND THE OTHER IN UPPER RGG 

TRIANGULAR FORM. IT REDUCES THE HESSENBERG MATRIX TORGG 
QUASI -TRIANGULAR FORM USING ORTHOGONAL TRANS FORMAT IONS RGG 
WHILE MAINTAINING THE TRIANGULAR FORM OF THE OTHER RGG 
MATRIX. RGG 

RGG 

QZVAL (PF2 62 ) ACCEPTS A PAIR OF REAL MATRICES, ONE OFRGG 

THEM IN QUASI -TRIANGULAR FORM AND THE OTHER IN UPPERRGG 

TRIANGULAR FORM. IT REDUCES THE QUASI -TRIANGULA RRGG 
MATRIX FURTHER, SO THAT ANY REMAINING 2-BY-2 BLOCKSRGG 
CORRESPOND TO PAIRS OF COMPLEX EIGENVALUES, AND RETURNSRGG 
QUANTITIES WHOSE RATIOS GIVE THE GENERALI 2 EDRGG 
EIGENVALUES. RGG 

RGG 

QZVEC ( PF2 63 ) ACCEPTS A PAIR OF REAL MATRICES, ONE OFRGG 
THEM IN QUASI -TRIANGULAR FORM (IN WHICH EACH 2-BY-2RGG 
BLOCK CORRESPONDS TO A PAIR OF COMPLEX EIGENVALUES) ANDRGG 
THE OTHER IN UPPER TRIANGULAR FORM. IT COMPUTES THERGG 
EIGENVECTORS OF THE TRIANGULAR PROBLEM AND TRANSFORMSRGG 
THE RESULTS BACK TO THE ORIGINAL COORDINATE SYSTEM. RGG 

RGG 
RGG 
RGG 
RGG 

PROGRAM TRGG (OUTPUT, TAPE6=OUTPUT) RGG 

RGG 

RGG 

DIMENSION A ( 5 , 5 ) ,B(5,5) , ALFR ( 5 ) ,ALFI(5) , BETA (5) ,Z(5,5) RGG 

RGG 


EXAMPLE 


N = 5 
NM = 5 
MAT2 


= 1 


DATA A /10. ,2. ,3. ,2*1. ,2, 
k 1. ,-1. ,1. ,2. ,1. ,9- 


, 12 . , 1 . , 2 . , 1 . , 3 , 

,3*1. ,-l. ,1. ,15. 


, 1 < 


,11m 

/ 


DATA B /12. ,1. , — 1 • ,2. ,2*1. ,14. ,1. , - 1 . ,1. ,-l. ,1. , 

* 16. ,-l. ,1. ,2. ,-l. ,-l. ,12. ,-l. ,3*1. ,-l. ,11. 

CALL RGG (NM, N, A, B , ALFR, ALFI , BETA, MATZ , Z, I ERR) 

WRITE (6,99) I ERR 

WRITE ( 6 , 100 ) ALFR, ALFI, BETA, ( ( Z ( I , J) , 1=1 , 5) ,J=1,5) 
FORMAT (1H 1,7 HI ERR = ,14) 


c 

* 

8H0ALFI 

= /1H , 5 ( G8 . 2 , 2X) / 

c 

* 

8H0BETA 

= /1H , 5 (G8 . 2 , 2X) / 

c 

★ 

5H0Z = 

/5(1H , 5 (G8 . 2 , 2X) / ) ) 

c 

STOP 


c 

END 



c 




c 

OUTPUT 

j 


c 




c 

I ERR = 

0 


c 

ALFR = 



c 

15. 

7.2 

16. 10. 

c 

ALFI = 



c 

0. 

0. 

0. 0. 


8.6 


RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 


91 

92 

93 

94 

95 

96 

97 

98 

99 
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123 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


BETA = 


9.9 

17. 

14 . 

11. 

Z = 
.26 

— . 59E-01 

.23 

-.30 

-.85 

.39 

1.0 

-.69 

1.0 

1.0 

.85 

.88 

-1.0 

.83 

-.39 

.72 

-.45 

-.84 

.65 

1.0 


13. 

- 1.0 

.26 

. 54E-01 
-.46 

-. 19E-01 


SUBROUTINE diverg(NM,N, A,B, ALFR, ALFI , BETA, MAT Z , Z, IERR) 
implicit real*8 (a-h,o-z) 

INTEGER N,NM, IERR, MAT Z „ /OT _ _ t . 

REAL* 8 A(NM,N) , B (NM,N) , ALFR (N) , ALFI (N) , BETA (N) , Z (NM, N) 

LOGICAL TF 

zero = 0.0e+00 

IF (N .LE. NM) GO TO 10 

IERR = 10 * N 

GO TO 50 


10 IF (MATZ .NE. 0) GO TO 20 
c FIND EIGENVALUES ONLY 

TF = .FALSE. 

CALL QZHES (NM, N , A, B, TF, Z) 

CALL QZIT(NM,N, A, B, zero ,TF,Z,IERR) 

CALL QZVAL(NM,N, A,B, ALFR, ALFI , BETA, TF, Z) 

GO TO 50 

c FIND BOTH EIGENVALUES AND EIGENVECTORS 

20 TF = .TRUE. 

CALL QZHES (NM, N, A, B,TF , Z) 

CALL QZIT(NM,N,A,B, zero ,TF,Z,IERR) 

CALL QZVAL(NM,N, A, B, ALFR, ALFI, BETA, TF,Z) 

IF (IERR .NE. 0) GO TO 50 

CALL QZVEC(NM,N, A, B, ALFR, ALFI, BETA, Z) 

50 RETURN 

C** THIS PROGRAM VALID ON FTN4 AND FTN5 ** 

END 

subroutine gotoer 

write (6, 10) . 

10 f ormat (' there is an error in calculating subroutine ) 

return 

end 

c 
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